Supplemental materials for: Rolek, B.W., McClure, CJW, Dunn, L., Curti, M., … Ridgway’s Hawk IPM and PVA

Contact information:

Metadata, data, and scripts used in analyses can be found at https://github.com/The-Peregrine-Fund/XXXXX.

The full workflow below is visible as a html website at: https://the-peregrine-fund.github.io/XXXXX/.

A permanent archive and DOI is available at: https://zenodo.org/doi/XXXXX


1. Code for IPM and PVA

1.1 Submodels without Integration

library('nimble')
library('parallel')
library ('coda')
load("/bsuscratch/brianrolek/riha_ipm/data.rdata")

#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
#   Observations (po) = y  
#     1 seen first-year (age=0, just before 1st b-day)
#     2 seen nonbreeder
#     3 seen breeder
#     4 not seen
#   States (ps)
#     1 alive first-year
#     2 alive nonbreeder
#     3 alive breeder
#     4 dead
#   Groups
#     1 wild-born
#     2 translocated and hacked
###################################################
# PARAMETERS
#   phiFY: survival probability first year 
#   phiA: survival probability nonbreeders
#   phiB: survival probability breeders
#   psiFYB: recruitment probability from first-year to breeder
#   psiAB: recruitment probability from nonbreeders to breeder
#   psiBA: recruitment probability from breeder to nonbreeders 
#   pFY: resight probability first-years
#   pA: resight probability nonbreeders
#   pB: resight probability breeders
#   lmu.prod: mean productivity (males and females) per territory on the log scale
#   sds: standard deviations for multivariate normal random effects for time and site
#   R: correlation coefficients for multivariate normal random effects for time and site
#   lambda: population growth rate (derived)
#   extinct: binary indicator of extirpation at a site
#   gamma: coefficient of effect from nest treatments
#   betas: coefficient of effect from translocations
#   deltas: coefficient of effect from survey effort
#   mus: overall means for survival, recruitment, and detections
#   r and rr: "r" parameter for negative binomial distribution
#             also described as omega in manuscript

#**********************
#* Model code
#**********************
mycode <- nimbleCode(
  {
    ###################################################
    # Priors and constraints
    ###################################################
    # Survival, recruitment, and detection can be correlated
    for (k in 1:8){ 
      betas[k] ~ dunif(-20, 20)  # prior for coefficients
    } # k
    for (j in 1:8){ 
      for (s in 1:nsite){
        lmus[j,s] <- logit(mus[j,s])
        mus[j,s] ~ dbeta(1,1) # prior for means
      }}     # m population #s sex #h hacked 
    
    # # Temporal random effects and correlations among all sites, synchrony
    # for (j in 1:p){ sds[j] ~ dexp(1) }# prior for temporal variation estimated using the multivariate normal distribution
    # R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
    # Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
    # U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p])
    # # Multivariate normal for temporal variance
    # for (t in 1:nyr){ 
    #   eps[1:p,t] ~ dmnorm(mu.zeroes[1:p],
    #                       cholesky = U[1:p, 1:p], prec_param = 0)
    # }
    
    # Temporal random effects and correlations between sites
    for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal variation estimated using the multivariate normal distribution
    R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
    Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
    U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2])
    # multivariate normal for temporal and site variance
    for (t in 1:nyr){ 
      for (s in 1:nsite){
        eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
                               cholesky = U2[1:p2, 1:p2], prec_param = 0)
      } } # s t 

    ###############################
    # Likelihood for productivity
    ###############################
    # Priors for productivity
    for (s in 1:nsite){
      lmu.f[s] ~ dnorm(0, sd=5)
    } # s
    gamma ~ dunif(-20, 20)
    rr ~ dexp(0.05)
    
    # Productivity likelihood      
    for (k in 1:npairsobs){
      f[k] ~ dnegbin(ppp[k], rr)
      ppp[k] <- rr/(rr+mu.f[k])
      log(mu.f[k]) <- lmu.f[site.pair[k]] +  
        gamma*treat.pair[k] + 
        #eps[9, year.pair[k] ] + 
        eta[9, site.pair[k], year.pair[k] ] 
    } # k
    # Derive yearly brood size for population model
    # Need to reorder because nimble doesn't 
    # handle nonconsecutive indices
    # yrind.pair is a matrix of indices for each site
    for (t in 1:nyr){
      for (s in 1:nsite){
        for (xxx in 1:pair.end[t,s]){
          fecmat[t,s,xxx] <- mu.f[ yrind.pair[xxx,t,s] ]
        } # xxx
        mn.f[t,s] <- mean( fecmat[t,s,1:pair.end[t,s]] )
      }} # s t
    
    # GOF for number of fledglings
    for (k in 1:npairsobs){
      f.obs[k] <- f[k] # observed counts
      f.exp[k] <- mu.f[k] # expected counts adult breeder
      f.rep[k] ~ dnegbin(ppp[k], rr) # expected counts
      f.dssm.obs[k] <- abs( ( f.obs[k] - f.exp[k] ) / (f.obs[k]+0.001) )
      f.dssm.rep[k] <- abs( ( f.rep[k] - f.exp[k] ) / (f.rep[k]+0.001) )
    } # k
    f.dmape.obs <- sum(f.dssm.obs[1:npairsobs])
    f.dmape.rep <- sum(f.dssm.rep[1:npairsobs])
    f.tvm.obs <- sd(brood[1:npairsobs])^2/mean(brood[1:npairsobs])
    f.tvm.rep <- sd(f.rep[1:npairsobs])^2/mean(f.rep[1:npairsobs])
    
    ################################
    # Likelihood for counts
    ################################
# Omitted. See IPM or PVA.
    
    ################################
    # Likelihood for survival
    ################################ 
    # Calculate yearly averages for sites for integration
    for (t in 1:nyr){
      for (s in 1:nsite){
        for (xxxx in 1:surv.end[t,s]){
          # Need to reorder because nimble doesn't 
          # handle nonconsecutive indices
          # yrind.surv is a matrix of indices for each site
          phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
          phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
          phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
          psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
          psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
          psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
          pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
          pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
        } # xxxx
        mn.phiFY[ t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] ) 
        mn.phiA[ t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
        mn.phiB[ t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiFYB[ t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiAB[ t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiBA[ t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
        mn.pA[ t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
        mn.pB[ t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
      }}
    
    for (i in 1:nind){
      for (t in 1:nyr){
        #Survival
        logit(phiFY[i,t]) <- eta[1, site[i,t],t] + # eps[1,t] + 
          lmus[1, site[i,t]] + betas[1]*hacked[i]  # first year
        logit(phiA[i,t]) <- eta[2, site[i,t],t] +# eps[2,t] + 
          lmus[2, site[i,t]] +  betas[2]*hacked[i] # nonbreeder
        logit(phiB[i,t]) <- eta[3, site[i,t],t] +# eps[3,t] + 
          lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
        #Recruitment
        logit(psiFYB[i,t]) <- eta[4, site[i,t],t] +# eps[4,t] + 
          lmus[4, site[i,t]] + betas[4]*hacked[i] # first year to breeder
        logit(psiAB[i,t]) <- eta[5, site[i,t],t] +# eps[5,t] + 
          lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
        logit(psiBA[i,t]) <- #eta[6, site[i,t],t] + eps[6,t] + 
          lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
        #Re-encounter
        logit(pA[i,t]) <- eta[7, site[i,t],t] +# eps[7,t] + 
          lmus[7, site[i,t]] + betas[7]*hacked[i] # resight of nonbreeders
        logit(pB[i,t]) <- eta[8, site[i,t],t] +# eps[8,t] + 
          lmus[8, site[i,t]] + betas[8]*hacked[i] # resight of breeders
      }#t
    }#i
    
    # Define state-transition and observation matrices
    for (i in 1:nind){  
      # Define probabilities of state S(t+1) given S(t)
      for (t in first[i]:(nyr-1)){
        ps[1,i,t,1] <- 0
        ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t]) 
        ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
        ps[1,i,t,4] <- (1-phiFY[i,t])
        
        ps[2,i,t,1] <- 0
        ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
        ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
        ps[2,i,t,4] <- (1-phiA[i,t])
        
        ps[3,i,t,1] <- 0
        ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
        ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
        ps[3,i,t,4] <- (1-phiB[i,t])
        
        ps[4,i,t,1] <- 0
        ps[4,i,t,2] <- 0
        ps[4,i,t,3] <- 0
        ps[4,i,t,4] <- 1
        
        # Define probabilities of O(t) given S(t)
        po[1,i,t,1] <- 1 
        po[1,i,t,2] <- 0
        po[1,i,t,3] <- 0
        po[1,i,t,4] <- 0
        
        po[2,i,t,1] <- 0
        po[2,i,t,2] <- pA[i,t]
        po[2,i,t,3] <- 0
        po[2,i,t,4] <- (1-pA[i,t])
        
        po[3,i,t,1] <- 0
        po[3,i,t,2] <- 0
        po[3,i,t,3] <- pB[i,t]
        po[3,i,t,4] <- (1-pB[i,t])
        
        po[4,i,t,1] <- 0
        po[4,i,t,2] <- 0
        po[4,i,t,3] <- 0
        po[4,i,t,4] <- 1
      } #t
    } #i
    
    # Likelihood 
    for (i in 1:nind){
      # Define latent state at first capture
      z[i,first[i]] <- y[i,first[i]]
      for (t in (first[i]+1):nyr){
        # State process: draw S(t) given S(t-1)
        z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
        # Observation process: draw O(t) given S(t)
        y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
      } #t
    } #i
  } )

#**********************
#* Function to run model in NIMBLE
#**********************
run_ipm <- function(info, datl, constl, code){
  library('nimble')
  library('coda')
  # helper function for multivariate normal
  uppertri_mult_diag <- nimbleFunction(
    run = function(mat = double(2), vec = double(1)) {
      returnType(double(2))
      p <- length(vec)
      out <- matrix(nrow = p, ncol = p, init = FALSE)
      for(i in 1:p){
        out[ ,i] <- mat[ ,i] * vec[i]
      }
      return(out)
    })
  assign('uppertri_mult_diag', uppertri_mult_diag, envir = .GlobalEnv)
  
  params <- c(# pop growth 
    #"lambda",
    # fecundity
    "lmu.f", "gamma", "rr", "mn.f", 
    # survival 
    "mus", "lmus", "betas",
    # abundance
    # "NB", "NF", "NFY", "N",
    # "r",
    # "N", "Ntot",
    # error terms
    #"eps", "sds", "Ustar", "U", "R",
    "eta", "sds2", "Ustar2", "U2", "R2",
    # yearly summaries
    'mn.phiFY', 'mn.phiA', 'mn.phiB',
    'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
    'mn.pA', 'mn.pB',
    # goodness of fit
    "f.dmape.obs", "f.dmape.rep",
    "f.tvm.obs", "f.tvm.rep"
    # "dmape.obs", "dmape.rep",
    # "tvm.obs", "tvm.rep",
    # "tturn.obs", "tturn.rep"
  )
  
  #n.chains=1; n.thin=200; n.iter=500000; n.burnin=300000
  #n.chains=1; n.thin=50; n.iter=100000; n.burnin=50000
  #n.chains=1; n.thin=10; n.iter=20000; n.burnin=10000
  n.chains=1; n.thin=200; n.iter=500000; n.burnin=300000
  
  mod <- nimbleModel(code, 
                     constants = constl, 
                     data = datl, 
                     inits = info$inits, 
                     buildDerivs = FALSE, # doesn't work when TRUE, no hope for HMC
                     calculate=T 
  ) 
  
  cmod <- compileNimble(mod, showCompilerOutput = TRUE)
  confhmc <- configureMCMC(mod)
  confhmc$setMonitors(params)
  hmc <- buildMCMC(confhmc)
  chmc <- compileNimble(hmc, project = mod, 
                        resetFunctions = TRUE,
                        showCompilerOutput = TRUE)
  
  post <- runMCMC(chmc,
                  niter = n.iter, 
                  nburnin = n.burnin,
                  nchains = n.chains,
                  thin = n.thin,
                  samplesAsCodaMCMC = T, 
                  setSeed = info$seed)
  
  return(post)
} # run_ipm function end

#*****************
#* Run chains in parallel
#*****************
this_cluster <- makeCluster(10)
post <- parLapply(cl = this_cluster, 
                  X = par_info, 
                  fun = run_ipm, 
                  datl = datl, 
                  constl = constl, 
                  code = mycode)
stopCluster(this_cluster)

save(post, mycode,
     file="/bsuscratch/brianrolek/riha_ipm/outputs/ipm_simp.rdata")

# save(post, mycode,
#      file="C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_statespace.rdata")

1.2 Integrated Population Model

library('nimble')
library('parallel')
library ('coda')
load("/bsuscratch/brianrolek/riha_ipm/data.rdata")

#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
#   Observations (po) = y  
#     1 seen first-year (age=0, just before 1st b-day)
#     2 seen nonbreeder
#     3 seen breeder
#     4 not seen
#   States (ps)
#     1 alive first-year
#     2 alive nonbreeder
#     3 alive breeder
#     4 dead
#   Groups
#     1 wild-born
#     2 translocated and hacked
###################################################
# PARAMETERS
#   phiFY: survival probability first year 
#   phiA: survival probability nonbreeders
#   phiB: survival probability breeders
#   psiFYB: recruitment probability from first-year to breeder
#   psiAB: recruitment probability from nonbreeders to breeder
#   psiBA: recruitment probability from breeder to nonbreeders 
#   pFY: resight probability first-years
#   pA: resight probability nonbreeders
#   pB: resight probability breeders
#   lmu.prod: mean productivity (males and females) per territory on the log scale
#   sds: standard deviations for multivariate normal random effects for time and site
#   R: correlation coefficients for multivariate normal random effects for time and site
#   lambda: population growth rate (derived)
#   extinct: binary indicator of extirpation at a site
#   gamma: coefficient of effect from nest treatments
#   betas: coefficient of effect from translocations
#   deltas: coefficient of effect from survey effort
#   mus: overall means for survival, recruitment, and detections
#   r and rr: "r" parameter for negative binomial distribution
#             also described as omega in manuscript

#**********************
#* Model code
#**********************
mycode <- nimbleCode(
  {
    ###################################################
    # Priors and constraints
    ###################################################
    # survival, recruitment, and detection can be correlated
    for (k in 1:8){
      betas[k] ~ dunif(-20, 20)  # prior for translocations coefficients
      deltas[k] ~ dunif(-20, 20) # prior for survey effort coefficients
    } # k
    
    for (j in 1:8){
      for (s in 1:nsite){
        lmus[j,s] <- logit(mus[j,s])
        mus[j,s] ~ dbeta(1,1) # prior for overall means
      }}     # 
    
    # # Temporal random effects and correlations among all sites
    # for (j in 1:p){ sds[j] ~ dexp(1) }# prior for temporal variation
    # # estimated using the multivariate normal distribution
    # R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
    # Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
    # U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p])
    # # multivariate normal for temporal variance
    # for (t in 1:nyr){ # survival params only have nyr-1, no problem to simulate from however
    #   eps[1:p,t] ~ dmnorm(mu.zeroes[1:p],
    #                       cholesky = U[1:p, 1:p], prec_param = 0)
    # }
    
    # Temporal random effects and correlations between sites
    for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal variation
    # estimated using the multivariate normal distribution
    R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
    Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
    U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2])
    # multivariate normal for temporal variance
    for (t in 1:nyr){ # survival params only have nyr-1, no problem to simulate from however
      for (s in 1:nsite){
        eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
                               cholesky = U2[1:p2, 1:p2], prec_param = 0)
      } } # s t 
    #######################
    # Derived params
    #######################
    for (s in 1:nsite){
      for (t in 1:(nyr-1)){
        lambda[t, s] <-  Ntot[t+1, s]/(Ntot[t, s])
        loglambda[t, s] <- log(lambda[t, s])
      }} #t
    
    ###############################
    # Likelihood for productivity
    ###############################
    # Priors
    for (s in 1:nsite){
      lmu.prod[s] ~ dnorm(0, sd=5)
    } # s
    gamma ~ dunif(-20, 20)
    rr ~ dexp(0.05)
    
    # Productivity likelihood      
    for (k in 1:npairsobs){
      prod[k] ~ dnegbin(ppp[k], rr)
      ppp[k] <- rr/(rr+mu.prod[k])
      log(mu.prod[k]) <- lmu.prod[site.pair[k]] +  
        gamma*treat.pair[k] + 
        #eps[9, year.pair[k] ] + 
        eta[9, site.pair[k], year.pair[k] ] 
    } # k
    # Derive yearly productivity for population model
    # need to reorder because nimble doesn't 
    # handle nonconsecutive indices
    # yrind.pair is a matrix of indices for each site
    for (t in 1:nyr){
      for (s in 1:nsite){
        for (xxx in 1:pair.end[t,s]){
          prodmat[t,s,xxx] <- mu.prod[ yrind.pair[xxx,t,s] ]
        } # xxx
        mn.prod[t,s] <- mean( prodmat[t,s,1:pair.end[t,s]] )
      }} # s t
    
    # GOF for productivity
    for (k in 1:npairsobs){
      f.obs[k] <- prod[k] # observed counts
      f.exp[k] <- mu.prod[k] # expected counts adult breeder
      f.rep[k] ~ dnegbin(ppp[k], rr) # expected counts
      f.dssm.obs[k] <- abs( ( f.obs[k] - f.exp[k] ) / (f.obs[k]+0.001) )
      f.dssm.rep[k] <- abs( ( f.rep[k] - f.exp[k] ) / (f.rep[k]+0.001) )
    } # k
    f.dmape.obs <- sum(f.dssm.obs[1:npairsobs])
    f.dmape.rep <- sum(f.dssm.rep[1:npairsobs])
    f.tvm.obs <- sd(brood[1:npairsobs])^2/mean(brood[1:npairsobs])
    f.tvm.rep <- sd(f.rep[1:npairsobs])^2/mean(f.rep[1:npairsobs])
    
    ################################
    # Likelihood for counts
    ################################
    # Abundance for year=1
    for (v in 1:7){ 
      for (s in 1:nsite){
        # subtract one to allow dcat to include zero
        N[v, 1, s] <- N2[v, 1, s] - 1 
        N2[v, 1, s] ~ dcat(pPrior[v, 1:s.end[v,s], s]) # Priors differ for FYs and adults
      }} # s t
    # Abundance for years > 1
    for (t in 1:(nyr-1)){
      for (s in 1:nsite){
        # Number of wild born juvs
        N[1, t+1, s] ~ dpois( (NFY[t, s]*mn.phiFY[t, s]*mn.psiFYB[t, s] + # first year breeders
                                 NF[t, s]*mn.phiA[t, s]*mn.psiAB[t, s] + # nonbreeders to breeders
                                 NB[t, s]*mn.phiB[t, s]*(1-mn.psiBA[t, s])) # breeders remaining
                              *mn.prod[t+1, s]/2 ) # end Poisson
        # Abundance of nonbreeders
        N[2, t+1, s] ~ dbin(mn.phiFY[t, s]*(1-mn.psiFYB[t, s]), NFY[t, s]) # Nestlings to nonbreeders
        N[3, t+1, s] ~ dbin(mn.phiA[t, s]*(1-mn.psiAB[t, s]), NF[t, s]) # Nonbreeders to nonbreeders
        N[4, t+1, s] ~ dbin(mn.phiB[t, s]*mn.psiBA[t, s], NB[t, s]) # Breeders to nonbreeders
        # Abundance of breeders
        N[5, t+1, s] ~ dbin(mn.phiFY[t, s]*mn.psiFYB[t, s], NFY[t, s]) # Nestlings to second year breeders
        N[6, t+1, s] ~ dbin(mn.phiA[t, s]*mn.psiAB[t, s], NF[t, s]) # Nonbreeder to breeder
        N[7, t+1, s] ~ dbin(mn.phiB[t, s]*(1-mn.psiBA[t, s]), NB[t, s]) # Breeder to breeder
      }} # s t
    
    # Count likelihoods, state-space model, and observation process    
    for (t in 1:nyr){
      for (s in 1:nsite){
        countsAdults[t, s] ~ dpois(lamAD[t, s]) # adult females 
        constraint_data[t, s] ~ dconstraint( (N[1, t, s] + hacked.counts[t, s]) >= 0 ) # constrain N1+hacked.counts to be >=0
        NFY[t, s] <- N[1, t, s] + hacked.counts[t, s] # Transfers translocated first-year females
        NF[t, s] <- sum(N[2:4, t, s]) # number of adult nonbreeder females
        NB[t, s] <- sum(N[5:7, t, s]) # number of adult breeder females
        NAD[t, s] <- NF[t, s] + NB[t, s] # number of adults
        Ntot[t, s] <- sum(N[1:7, t, s]) # total number of females
        log(lamAD[t, s]) <- log(NAD[t, s]) + deltas[1]*effort2[t, s] + deltas[2]*effort2[t, s]^2
        log(lamFY[t, s]) <- log(N[1, t, s]) + deltas[3]*effort2[t, s] + deltas[4]*effort2[t, s]^2
        }# s
      # First-years at different sites have different distributions
      # for better model fit
      countsFY[t, 1] ~ dpois(lamFY[t, 1]) # doesn't have any zeroes so poisson
      countsFY[t, 2] ~ dnegbin(pp[t], r) # first year females, includes translocated/hacked
      pp[t] <- r/(r+(lamFY[t, 2] ))
    } # t
    r ~ dexp(0.05)
    ###################
    # Assess GOF of the state-space models for counts
    # Step 1: Compute statistic for observed data
    # Step 2: Use discrepancy measure: mean absolute error
    # Step 3: Use test statistic: number of turns
    ###################
    for (t in 1:nyr){
      c.repFY[t, 1] ~ dpois( lamFY[t, 1] )
      c.repFY[t, 2] ~ dnegbin( pp[t], r )
      for (s in 1:nsite){
        c.expAD[t, s] <- lamAD[t, s]  # expected counts adult breeder
        c.expFY[t, s] <- lamFY[t, s]
        c.obsAD[t, s] <- countsAdults[t, s]
        c.obsFY[t, s] <- countsFY[t, s]  # first year
        c.repAD[t, s] ~ dpois( lamAD[t, s] ) # simulated counts
        dssm.obsAD[t, s] <- abs( ( (c.obsAD[t, s]) - (c.expAD[t, s]) ) / (c.obsAD[t, s]+0.001)  )
        dssm.obsFY[t, s] <- abs( ( (c.obsFY[t, s]) - (c.expFY[t, s]) ) / (c.obsFY[t, s]+0.001)  )
        dssm.repAD[t, s] <- abs( ( (c.repAD[t, s]) - (c.expAD[t, s]) ) / (c.repAD[t, s]+0.001) )
        dssm.repFY[t, s] <- abs( ( (c.repFY[t, s]) - (c.expFY[t, s]) ) / (c.repFY[t, s]+0.001) )
      }} # t
    dmape.obs[1] <- sum(dssm.obsAD[1:nyr, 1:nsite])
    dmape.obs[2] <- sum(dssm.obsFY[1:nyr, 1:nsite])
    dmape.rep[1] <- sum(dssm.repAD[1:nyr, 1:nsite])
    dmape.rep[2] <- sum(dssm.repFY[1:nyr, 1:nsite])
    tvm.obs[1] <- sd(c.obsB[1:nyr, 1:nsite])^2/mean(c.obsB[1:nyr, 1:nsite])
    tvm.obs[2] <- sd(c.obsFY[1:nyr, 1:nsite])^2/mean(c.obsFY[1:nyr, 1:nsite])
    tvm.rep[1] <- sd(c.repB[1:nyr, 1:nsite])^2/mean(c.repB[1:nyr, 1:nsite])
    tvm.rep[2] <- sd(c.repFY[1:nyr, 1:nsite])^2/mean(c.repFY[1:nyr, 1:nsite])
    # # Test statistic for number of turns
    for (s in 1:nsite){
      for (t in 1:(nyr-2)){
        tt1.obsAD[t, s] <- step(c.obsAD[t+2, s] - c.obsAD[t+1, s])
        tt2.obsAD[t, s] <- step(c.obsAD[t+1, s] - c.obsAD[t, s])
        tt3.obsAD[t, s] <- equals(tt1.obsAD[t, s] + tt2.obsAD[t, s], 1)
        tt1.obsFY[t, s] <- step(c.obsFY[t+2, s] - c.obsFY[t+1, s])
        tt2.obsFY[t, s] <- step(c.obsFY[t+1, s] - c.obsFY[t, s])
        tt3.obsFY[t, s] <- equals(tt1.obsFY[t, s] + tt2.obsFY[t, s], 1)
      }} # t
    tturn.obs[1] <- sum(tt3.obsAD[1:(nyr-2), 1:nsite])
    tturn.obs[2] <- sum(tt3.obsFY[1:(nyr-2), 1:nsite])
    for (s in 1:nsite){
      for (t in 1:(nyr-2)){
        tt1.repAD[t, s] <- step(c.repAD[t+2, s] - c.repAD[t+1, s])
        tt2.repAD[t, s] <- step(c.repAD[t+1, s] - c.repAD[t, s])
        tt3.repAD[t, s] <- equals(tt1.repAD[t, s] + tt2.repAD[t, s], 1)
        tt1.repFY[t, s] <- step(c.repFY[t+2, s] - c.repFY[t+1, s])
        tt2.repFY[t, s] <- step(c.repFY[t+1, s] - c.repFY[t, s])
        tt3.repFY[t, s] <- equals(tt1.repFY[t, s] + tt2.repFY[t, s], 1)
      }} # t
    tturn.rep[1] <- sum(tt3.repAD[1:(nyr-2), 1:nsite])
    tturn.rep[2] <- sum(tt3.repFY[1:(nyr-2), 1:nsite])
    
    ################################
    # Likelihood for survival
    ################################ 
    # Calculate averages for sites each year for integration
    for (t in 1:nyr){
      for (s in 1:nsite){
        for (xxxx in 1:surv.end[t,s]){
          # Reorder because nimble doesn't 
          # handle nonconsecutive indices
          # yrind.surv is a matrix of indices for each site
          phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
          phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
          phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
          psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
          psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
          psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
          pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
          pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
        } # xxxx
        mn.phiFY[ t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] ) 
        mn.phiA[ t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
        mn.phiB[ t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiFYB[ t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiAB[ t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiBA[ t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
        mn.pA[ t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
        mn.pB[ t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
      }}
    
    for (i in 1:nind){
      for (t in 1:nyr){
        #Survival
        logit(phiFY[i,t]) <- eta[1, site[i,t],t] + # eps[1,t] + 
          lmus[1, site[i,t]] + betas[1]*hacked[i]  # first year
        logit(phiA[i,t]) <- eta[2, site[i,t],t] + #eps[2,t] + 
          lmus[2, site[i,t]] +  betas[2]*hacked[i] # nonbreeder
        logit(phiB[i,t]) <- eta[3, site[i,t],t] + #eps[3,t] + 
          lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
        #Recruitment
        logit(psiFYB[i,t]) <- eta[4, site[i,t],t] + #eps[4,t] + 
          lmus[4, site[i,t]] + betas[4]*hacked[i] # first year to breeder
        logit(psiAB[i,t]) <- eta[5, site[i,t],t] + #eps[5,t] + 
          lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
        logit(psiBA[i,t]) <- #eta[6, site[i,t],t] + eps[6,t] + 
          lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
        #Re-encounter
        logit(pA[i,t]) <- eta[7, site[i,t],t] + #eps[7,t] + 
          lmus[7, site[i,t]] + betas[7]*hacked[i] +
          deltas[5]*effort2[t, site[i,t]] + deltas[6]*effort2[t, site[i,t]]^2# resight of nonbreeders
        logit(pB[i,t]) <- eta[8, site[i,t],t] + #eps[8,t] + 
          lmus[8, site[i,t]] + betas[8]*hacked[i] + 
          deltas[7]*effort2[t, site[i,t]] + deltas[8]*effort2[t, site[i,t]]^2# resight of breeders
      }#t
    }#i
    
    # Define state-transition and observation matrices
    for (i in 1:nind){  
      # Define probabilities of state S(t+1) given S(t)
      for (t in first[i]:(nyr-1)){
        ps[1,i,t,1] <- 0
        ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t]) 
        ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
        ps[1,i,t,4] <- (1-phiFY[i,t])
        
        ps[2,i,t,1] <- 0
        ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
        ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
        ps[2,i,t,4] <- (1-phiA[i,t])
        
        ps[3,i,t,1] <- 0
        ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
        ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
        ps[3,i,t,4] <- (1-phiB[i,t])
        
        ps[4,i,t,1] <- 0
        ps[4,i,t,2] <- 0
        ps[4,i,t,3] <- 0
        ps[4,i,t,4] <- 1
        
        # Define probabilities of O(t) given S(t)
        po[1,i,t,1] <- 1 
        po[1,i,t,2] <- 0
        po[1,i,t,3] <- 0
        po[1,i,t,4] <- 0
        
        po[2,i,t,1] <- 0
        po[2,i,t,2] <- pA[i,t]
        po[2,i,t,3] <- 0
        po[2,i,t,4] <- (1-pA[i,t])
        
        po[3,i,t,1] <- 0
        po[3,i,t,2] <- 0
        po[3,i,t,3] <- pB[i,t]
        po[3,i,t,4] <- (1-pB[i,t])
        
        po[4,i,t,1] <- 0
        po[4,i,t,2] <- 0
        po[4,i,t,3] <- 0
        po[4,i,t,4] <- 1
      } #t
    } #i
    
    # Likelihood 
    for (i in 1:nind){
      # Define latent state at first capture
      z[i,first[i]] <- y[i,first[i]]
      for (t in (first[i]+1):nyr){
        # State process: draw S(t) given S(t-1)
        z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
        # Observation process: draw O(t) given S(t)
        y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
      } #t
    } #i
  } )

#**********************
#* Function to run model in NIMBLE
#**********************
run_ipm <- function(info, datl, constl, code){
  library('nimble')
  library('coda')
  # helper function for multivariate normal
  uppertri_mult_diag <- nimbleFunction(
    run = function(mat = double(2), vec = double(1)) {
      returnType(double(2))
      p <- length(vec)
      out <- matrix(nrow = p, ncol = p, init = FALSE)
      for(i in 1:p){
        out[ ,i] <- mat[ ,i] * vec[i]
      }
      return(out)
    })
  assign('uppertri_mult_diag', uppertri_mult_diag, envir = .GlobalEnv)
  
  params <- c(
    # pop growth 
    "lambda",
    # fecundity
    "lmu.prod", "gamma", "rr", "mn.prod", 
    # survival 
    "mus", "lmus", "betas", "deltas",
    # abundance
    "NB", "NF", "NFY", "N", "NAD",
    "r",
    "N", "Ntot",
    # error terms
    #"eps", "sds", "Ustar", "U", "R",
    "eta", "sds2", "Ustar2", "U2", "R2",
    # yearly summaries
    'mn.phiFY', 'mn.phiA', 'mn.phiB',
    'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
    'mn.pA', 'mn.pB',
    # goodness of fit
    "f.dmape.obs", "f.dmape.rep",
    "f.tvm.obs", "f.tvm.rep",
    "dmape.obs", "dmape.rep",
    "tvm.obs", "tvm.rep",
    "tturn.obs", "tturn.rep"
  )
  #n.chains=1; n.thin=1; n.iter=500; n.burnin=100
  #n.chains=1; n.thin=200; n.iter=800000; n.burnin=600000
  n.chains=1; n.thin=200; n.iter=600000; n.burnin=400000
  
  mod <- nimbleModel(code, 
                     constants = constl, 
                     data = datl, 
                     inits = info$inits, 
                     buildDerivs = FALSE, # doesn't work when TRUE, no hope for HMC
                     calculate=T 
  ) 
  
  cmod <- compileNimble(mod, showCompilerOutput = TRUE)
  confhmc <- configureMCMC(mod)
  confhmc$setMonitors(params)
  hmc <- buildMCMC(confhmc)
  chmc <- compileNimble(hmc, project = mod, 
                        resetFunctions = TRUE,
                        showCompilerOutput = TRUE)
  
  post <- runMCMC(chmc,
                  niter = n.iter, 
                  nburnin = n.burnin,
                  nchains = n.chains,
                  thin = n.thin,
                  samplesAsCodaMCMC = T, 
                  setSeed = info$seed)
  
  return(post)
} # run_ipm function end

#*****************
#* Run chains in parallel
#*****************
par_info <- par_info[1:10]
this_cluster <- makeCluster(10)
post <- parLapply(cl = this_cluster, 
                  X = par_info, 
                  fun = run_ipm, 
                  datl = datl, 
                  constl = constl, 
                  code = mycode)
stopCluster(this_cluster)

# save(post, mycode,
#      file="/bsuscratch/brianrolek/riha_ipm/outputs/ipm_shortrun.rdata")

save(post, mycode,
     file="/bsuscratch/brianrolek/riha_ipm/outputs/ipm_longrun.rdata")

1.3 Population Viability Analysis

library('nimble')
library('parallel')
library ('coda')

load("/bsuscratch/brianrolek/riha_ipm/data.rdata")
source("/bsuscratch/brianrolek/riha_ipm/MCMCvis.R")
load("/bsuscratch/brianrolek/riha_ipm/outputs/ipm_longrun.rdata")

# load("C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//outputs//ipm_longrun.rdata")
# load("data//data.rdata")
# library ("MCMCvis")

#**********************
#* Set data to reflect PVA scenarios
#**********************
constl$K <- 50 # number of future years
cpus <- 10 # number of processors
constl$SC <- 36 # number of PVA scenarios
datl$constraint_data <- rbind(datl$constraint_data, array(1, dim=c(constl$K,2)) ) # to help constrain FYs born + hacked to be positive
constl$effort2 <- rbind(constl$effort2, array(0, dim=c(constl$K,2), dimnames=list(2024:(2024+constl$K-1), c("LH", "PC"))))

constl$num.treated <- rep( c(0, 15, 30, 45, 0, 15, 30, 45, 0, 15, 30, 45), each=3 ) # 100s sub in for All but are over-ridden in model code
num.hacked <- rep( c(0, 0, 0, 0, 5, 5, 5, 5, 10, 10, 10, 10), each=3 )
surv.diff <- array(NA, dim=c(constl$SC, 3, 2), dimnames=list(scenario=1:constl$SC, stage=c('FY', 'NB', 'B'), site=c('LH', 'PC')))
surv.diff[,1,] <- matrix(c(rep( c(0, 0.130, 0.260), 12),
                          rep( c(0, 0.146, 0.292), 12)), ncol=2)
surv.diff[,2,] <- matrix(c(rep( c(0, 0.006, 0.012), 12),
                           rep( c(0, 0.022, 0.044), 12)), ncol=2)
surv.diff[,3,] <- matrix(c(rep( c(0, 0.016, 0.032), 12),
                           rep( c(0, 0.026, 0.052), 12)), ncol=2)
constl$surv.diff <- surv.diff 

hacked.counts <- array(0, dim=c(constl$SC, constl$nyr+constl$K, 2))
for (sc in 1:constl$SC){
  hacked.counts[sc,1:constl$nyr, 1:2] <- constl$hacked.counts[, 1:2] 
  hacked.counts[sc,14:(constl$nyr+constl$K), 1] <- -num.hacked[sc]
  hacked.counts[sc,14:(constl$nyr+constl$K), 2] <- 0
}
constl$hacked.counts <- hacked.counts

#**********************
#* Specify priors taken from posterior of IPM
#**********************
#* This helps ensure that all chains run
# Identify chains with NAs that failed to initialize
out <- lapply(post, as.mcmc)
NAlist <- c()
for (i in 1:length(out)){
  NAlist[i] <- any (is.na(out[[i]][,1:286]) | out[[i]][,1:286]<0)
}
out <- out[!NAlist]
outp <- MCMCpstr(out, type="chains")
!NAlist

Ni.func <- function (){
  # randomly select inits from ones that work
  # Ni <- outp$N[1:7,1:constl$nyr,1:2,
  #              sample(seq(1, 10000, by=100), 1, replace = F)]
  # take the means of the posterior for inits
  Ni <- apply(outp$N, c(1,2,3), mean) |> round()
  Ni.pva <- array(NA, dim=c(constl$SC, dim(Ni)+c(0,constl$K,0)))
  for (sc in 1:constl$SC){
    Ni.pva[sc, 1:7, 1:constl$nyr, 1:2] <- Ni
    for (t in constl$nyr:(constl$nyr+constl$K)){
      for (s in 1:2){
      Ni.pva[sc, 1:7, t, s] <- Ni[1:7, 13, s] 
    }}} # t sc
  return(Ni.pva)
} # function

u2 <- apply(outp$Ustar2, c(1,2), mean)

inits.func.pva <- function (){
  list(  
    # fecundity 
    lmu.prod = apply(outp$lmu.prod, 1, mean),
    gamma = mean(outp$gamma), 
    rr = mean(outp$rr),
    # survival
    z = par_info[[1]]$inits$z, 
    mus = apply(outp$mus, c(1,2), mean), 
    betas = apply(outp$betas, 1, mean),
    deltas = apply(outp$deltas, 1, mean),
    sds2 =  apply(outp$sds2, 1, mean),
    Ustar2 = u2,
    # counts
    countsAdults= matrix(c(374, 335, 305, 295, rep(NA, length(2015:2023)), rep(NA, length(2011:2023)) ), nrow=13), 
    r = mean(outp$r),
    N = Ni.func()
  )}

# Set seed then save for reproducibility
# then randomly draw for each chain
set.seed(1)
seeds <- sample(1:1000000, size=cpus, replace=FALSE)
par_info_pva <- list()
for (i in 1:cpus){
  par_info_pva[[i]] <- list(seed=seeds[i], inits = inits.func.pva())
}

#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
#   Observations (po) = y  
#     1 seen first-year (age=0, just before 1st b-day)
#     2 seen nonbreeder
#     3 seen breeder
#     4 not seen
#   States (ps)
#     1 alive first-year
#     2 alive nonbreeder
#     3 alive breeder
#     4 dead
#   Groups
#     1 wild-born
#     2 translocated and hacked
###################################################
# PARAMETERS
#   phiFY: survival probability first year 
#   phiA: survival probability nonbreeders
#   phiB: survival probability breeders
#   psiFYB: recruitment probability from first-year to breeder
#   psiAB: recruitment probability from nonbreeders to breeder
#   psiBA: recruitment probability from breeder to nonbreeders 
#   pFY: resight probability first-years
#   pA: resight probability nonbreeders
#   pB: resight probability breeders
#   lmu.prod: mean productivity (males and females) per territory on the log scale
#   sds: standard deviations for multivariate normal random effects for time and site
#   R: correlation coefficients for multivariate normal random effects for time and site
#   lambda: population growth rate (derived)
#   extinct: binary indicator of extirpation at a site
#   gamma: coefficient of effect from nest treatments
#   betas: coefficient of effect from translocations
#   deltas: coefficient of effect from survey effort
#   mus: overall means for survival, recruitment, and detections
#   r and rr: "r" parameter for negative binomial distribution
#             also described as omega in manuscript

#**********************
#* Model code
#**********************
mycode <- nimbleCode(
  {
    ###################################################
    # Priors and constraints
    ###################################################
    # survival, recruitment, and detection can be correlated
    for (k in 1:8){ 
      betas[k] ~ dunif(-20, 20)  # prior for coefficients
      deltas[k] ~ dunif(-20, 20)
    } # k
    
    for (j in 1:8){
      for (s in 1:nsite){
        lmus[j,s] <- logit(mus[j,s])
        mus[j,s] ~ dbeta(1,1) # prior for means
      }}     # m population #s sex #h hacked
    
    # Temporal and site random effects using the multivariate normal distribution
    for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal and site variation
    R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
    Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
    U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2]) # custom function taken from nimble listserv (specified below)
    # multivariate normal for temporal variance
    for (t in 1:(nyr+K) ){ 
      for (s in 1:nsite){
        eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
                               cholesky = U2[1:p2, 1:p2], prec_param = 0)
      } } # s t 
    #######################
    # Derived params
    #######################
    for(sc in 1:SC){
    for (s in 1:nsite){
      for (t in 1:(nyr-1+K)){
        lambda[sc, t, s] <-  Ntot[sc, t+1, s]/(Ntot[sc, t, s]) # population groewth rate
        loglambda[sc, t, s] <- log(lambda[sc, t, s]) # r
      }} #t
    
    for (s in 1:nsite){
      for (t in 1:K){
        # quasi-extinction probabilities given population segments, N total, adults, and breeders
        extinct[sc, t, s] <- equals(Ntot[sc, nyr+t, s], 0) 
        extinctAD[sc, t, s] <- equals(NAD[sc, nyr+t, s], 0)
        extinctB[sc, t, s] <- equals(NB[sc, nyr+t, s], 0)
      }}
    } # sc
    
    ###############################
    # Likelihood for productivity
    ###############################
    # Priors 
    for (s in 1:nsite){
      lmu.prod[s] ~ dnorm(0, sd=5)
    } # s
    gamma ~ dunif(-20, 20)
    rr ~ dexp(0.05)
    
    # Productivity likelihood       
    for (k in 1:npairsobs){
      prod[k] ~ dnegbin(ppp[k], rr)
      ppp[k] <- rr/(rr+mu.prod[k])
      log(mu.prod[k]) <- lmu.prod[site.pair[k]] +  
        gamma*treat.pair[k] + 
        eta[9, site.pair[k], year.pair[k] ] 
    } # k
    # Derive yearly productivity for population model
    # need to reorder because nimble doesn't 
    # handle nonconsecutive indices
    # yrind.pair is a matrix of indices for each site
    for (t in 1:nyr){
      for (s in 1:nsite){
        for (xxx in 1:pair.end[t,s]){
          prodmat[t,s,xxx] <- mu.prod[ yrind.pair[xxx,t,s] ]
        } # xxx
        mn.prod[1,t,s] <- mean( prodmat[t,s,1:pair.end[t,s]] )
        for(sc in 2:SC){
          mn.prod[sc, t, s] <- mn.prod[1,t,s]
      }}} # s t
    # Future projections- fecundity
      for(sc in 1:SC){
        for (t in (nyr+1):(nyr+K) ){
          for (s in 1:nsite){
            perc.treat[sc, t, s] <-  # 
              step( NB[sc, t, s]-num.treated[sc] ) * (num.treated[sc] + 0.001)/(NB[sc, t, s]+0.001)  + # num.treated > NB, calculate proportion
              1-step( NB[sc, t, s]-num.treated[sc] )   #  num.treated < NB set to 1
            
        log(mn.prod[sc, t, s]) <- lmu.prod[s] + 
          gamma*perc.treat[sc, t, s] + # treat.pair set to one here.
                              eta[9, s, t] 
  }} } # s t sc
    
    ################################
    # Likelihood for counts
    ################################
    # Abundance for year=1
      for (v in 1:7){ 
      for (s in 1:nsite){
        # Subtract one to allow dcat to include zero
        N[1, v, 1, s] <- N2[1, v, 1, s] - 1 
        N2[1, v, 1, s] ~ dcat(pPrior[v, 1:s.end[v,s], s])
        # Assign estimates for years with data to all scenarios
        for (sc in 2:SC){
          N[sc, v, 1, s] <- N[1, v, 1, s]
        }
      }} # s t

    # Abundance for years > 1
    for (t in 1:(nyr-1)){
      for (s in 1:nsite){
        # Number of wild born juvs
        N[1, 1, t+1, s] ~ dpois( (NFY[1, t, s]*mn.phiFY[1, t, s]*mn.psiFYB[1, t, s] + # first year breeders
                                 NF[1, t, s]*mn.phiA[1, t, s]*mn.psiAB[1, t, s] + # nonbreeders to breeders
                                 NB[1, t, s]*mn.phiB[1, t, s]*(1-mn.psiBA[1, t, s])) # breeders remaining
                              *(mn.prod[1, t+1, s]/2) ) # end Poisson
        # Abundance of nonbreeders
        N[1, 2, t+1, s] ~ dbin(mn.phiFY[1, t, s]*(1-mn.psiFYB[1, t, s]), NFY[1, t, s]) # Nestlings to second year nonbreeders
        N[1, 3, t+1, s] ~ dbin(mn.phiA[1, t, s]*(1-mn.psiAB[1, t, s]), NF[1, t, s]) # Nonbreeders to nonbreeders
        N[1, 4, t+1, s] ~ dbin(mn.phiB[1, t, s]*mn.psiBA[1, t, s], NB[1, t, s]) # Breeders to nonbreeders
        # Abundance of breeders
        N[1, 5, t+1, s] ~ dbin(mn.phiFY[1, t, s]*mn.psiFYB[1, t, s], NFY[1, t, s]) # Nestlings to second year breeders
        N[1, 6, t+1, s] ~ dbin(mn.phiA[1, t, s]*mn.psiAB[1, t, s], NF[1, t, s]) # Nonbreeder to breeder
        N[1, 7, t+1, s] ~ dbin(mn.phiB[1, t, s]*(1-mn.psiBA[1, t, s]), NB[1, t, s]) # Breeder to breeder
        # Assign estimates for years with data to 
        # all scenarios
        for (v in 1:7){ 
        for (sc in 2:SC){
          N[sc, v, t+1, s] <- N[1, v, t+1, s]
        } } # sc
        }} # s t

  # Future projections- population model 
    for (sc in 1:SC){
      for (t in nyr:(nyr+K-1)){
        for (s in 1:nsite){
          # Number of wild born juvs
          N[sc, 1, t+1, s] ~ dpois( (NFY[sc, t, s]*mn.phiFY[sc, t, s]*mn.psiFYB[sc, t, s] + # first year breeders
                                       NF[sc, t, s]*mn.phiA[sc, t, s]*mn.psiAB[sc, t, s] + # nonbreeders to breeders
                                       NB[sc, t, s]*mn.phiB[sc, t, s]*(1-mn.psiBA[sc, t, s])) # breeders remaining
                                    *(mn.prod[sc, t+1, s]/2) ) # end Poisson
          # Abundance of nonbreeders
          ## Second year nonbreeders
          N[sc, 2, t+1, s] ~ dbin(mn.phiFY[sc, t, s]*(1-mn.psiFYB[sc, t, s]), NFY[sc, t, s]) # Nestlings to second year nonbreeders
          ## Adult nonbreeders
          N[sc, 3, t+1, s] ~ dbin(mn.phiA[sc, t, s]*(1-mn.psiAB[sc, t, s]), NF[sc, t, s]) # Nonbreeders to nonbreeders
          N[sc, 4, t+1, s] ~ dbin(mn.phiB[sc, t, s]*mn.psiBA[sc, t, s], NB[sc, t, s]) # Breeders to nonbreeders
          # Abundance of breeders
          ## Second year breeders
          N[sc, 5, t+1, s] ~ dbin(mn.phiFY[sc, t, s]*mn.psiFYB[sc, t, s], NFY[sc, t, s]) # Nestlings to second year breeders
          ## Adult breeders
          N[sc, 6, t+1, s] ~ dbin(mn.phiA[sc, t, s]*mn.psiAB[sc, t, s], NF[sc, t, s]) # Nonbreeder to breeder
          N[sc, 7, t+1, s] ~ dbin(mn.phiB[sc, t, s]*(1-mn.psiBA[sc, t, s]), NB[sc, t, s]) # Breeder to breeder
        }} # s t
    } # sc
    
    # Count likelihoods, state-space model, and observation process 
    for (t in 1:nyr){
      for (s in 1:nsite){
        constraint_data[t, s] ~ dconstraint( (N[1, 1, t, s] + hacked.counts[1, t, s]) >= 0 ) # Transfers translocated first-year females
        NFY[1, t, s] <- N[1, 1, t, s] + hacked.counts[1, t, s] # Transfers translocated first-year females
        NF[1, t, s] <- sum(N[1, 2:4, t, s]) # number of adult nonbreeder females
        NB[1, t, s] <- sum(N[1, 5:7, t, s]) # number of adult breeder females
        NAD[1, t, s] <- NF[1, t, s] + NB[1, t, s] # number of adults
        Ntot[1, t, s] <- sum(N[1, 1:7, t, s]) # total number of females
        countsAdults[t, s] ~ dpois(lamAD[1, t, s]) # adult females
        # constrain N1+hacked.counts to be >=0
        log(lamAD[1, t, s]) <- log(NAD[1, t, s]) + deltas[1]*effort2[t, s] + deltas[2]*effort2[t, s]^2
        log(lamFY[1, t, s]) <- log(N[1, 1, t, s]) + deltas[3]*effort2[t, s] + deltas[4]*effort2[t, s]^2
        for (sc in 2:SC){
          lamAD[sc, t, s] <- lamAD[1, t, s]
          lamFY[sc, t, s] <- lamFY[1, t, s]
          NFY[sc, t, s] <- NFY[1, t, s]
          NF[sc, t, s] <- NF[1, t, s]
          NB[sc, t, s] <- NB[1, t, s]
          NAD[sc, t, s] <- NAD[1, t, s]
          Ntot[sc, t, s] <- Ntot[1, t, s]
        } # sc
      }# s
      # First-years at different sites have different distributions
      # for better model fit
      countsFY[t, 1] ~ dpois(lamFY[1, t, 1]) # doesn't have any zeroes so poisson
      countsFY[t, 2] ~ dnegbin(pp[t], r) # first year females, includes translocated/hacked
      pp[t] <- r/(r+(lamFY[1, t, 2] ))
    } # t
    r ~ dexp(0.05)
    
   # Future projections of counts 
    for(sc in 1:SC){
      for (t in (nyr+1):(nyr+K)){
        for (s in 1:nsite){
          # constrain N1+hacked.counts to be >=0, and allow it to shrink as the population shrinks
          NFY[sc, t, s] <- step( (N[sc, 1, t, s] + hacked.counts[sc, t, s]) ) * (N[sc, 1, t, s] + hacked.counts[sc, t, s])  
          NF[sc, t, s] <- sum(N[sc, 2:4, t, s]) # number of adult nonbreeder females
          NB[sc, t, s] <- sum(N[sc, 5:7, t, s]) # number of adult breeder females
          NAD[sc, t, s] <- NF[sc, t, s] + NB[sc, t, s] # number of adults
          Ntot[sc, t, s] <- sum(N[sc, 1:7, t, s]) # total number of females
        }# s
      }  }# t sc
    
    ################################
    # Likelihood for survival
    ################################ 
    # Calculate yearly averages for integration
    for (t in 1:(nyr-1)){
      for (s in 1:nsite){
        for (xxxx in 1:surv.end[t,s]){
          # need to reorder because nimble doesn't 
          # handle nonconsecutive indices
          # yrind.surv is a matrix of indices for each site
          phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
          phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
          phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
          psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
          psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
          psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
          pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
          pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
        } # xxxx
        mn.phiFY[1, t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] ) 
        mn.phiA[1, t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
        mn.phiB[1, t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiFYB[1, t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiAB[1, t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiBA[1, t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
        mn.pA[1, t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
        mn.pB[1, t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
        for (sc in 2:SC){
          mn.phiFY[sc, t, s] <- mn.phiFY[1, t, s] 
          mn.phiA[sc, t, s] <- mn.phiA[1, t, s] 
          mn.phiB[sc, t, s] <- mn.phiB[1, t, s] 
          mn.psiFYB[sc, t, s] <- mn.psiFYB[1, t, s]
          mn.psiAB[sc, t, s] <- mn.psiAB[1, t, s]
          mn.psiBA[sc, t, s] <- mn.psiBA[1, t, s]
          mn.pA[sc, t, s] <- mn.pA[1, t, s]
          mn.pB[sc, t, s] <- mn.pB[1, t, s]
        }
      }}
    
    for (i in 1:nind){
      for (t in 1:(nyr-1)){
        #Survival
        logit(phiFY[i,t]) <- eta[1, site[i,t],t] + 
          lmus[1, site[i,t]] + betas[1]*hacked[i]  # first year
        logit(phiA[i,t]) <- eta[2, site[i,t],t] + 
          lmus[2, site[i,t]] #+  betas[2]*hacked[i] # nonbreeder
        logit(phiB[i,t]) <- eta[3, site[i,t],t] +  
          lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
        #Recruitment
        logit(psiFYB[i,t]) <- eta[4, site[i,t],t] +  
          lmus[4, site[i,t]] #+ betas[4]*hacked[i] # first year to breeder
        logit(psiAB[i,t]) <- eta[5, site[i,t],t] + 
          lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
        logit(psiBA[i,t]) <-  
          lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
        #Re-encounter
        logit(pA[i,t]) <- eta[7, site[i,t],t] + 
          lmus[7, site[i,t]] + betas[7]*hacked[i] +
          deltas[5]*effort2[t, site[i,t]] + deltas[6]*effort2[t, site[i,t]]^2 # resight of nonbreeders
        logit(pB[i,t]) <- eta[8, site[i,t],t] + 
          lmus[8, site[i,t]] + #betas[8]*hacked[i] + 
          deltas[7]*effort2[t, site[i,t]] + deltas[8]*effort2[t, site[i,t]]^2 # resight of breeders
      }#t
    }#i
    
    # Future projections of survival, recruitment, and detection
    cap <- 0.97
    for(sc in 1:SC){ 
      for (t in nyr:(nyr+K)){
        for (s in 1:nsite){
        #Survival
        mn.phiFY1[sc, t, s] <- ilogit( eta[1, s, t] + 
          lmus[1, s] ) + surv.diff[sc, 1, s]  
        mn.phiA1[sc, t, s] <- ilogit( eta[2, s, t] + 
          lmus[2, s] ) + surv.diff[sc, 2, s]
        mn.phiB1[sc, t, s] <- ilogit( eta[3, s, t] +  
          lmus[3, s] ) + surv.diff[sc, 3, s]
        # enforce the cap
        mn.phiFY[sc, t, s] <- step(cap - mn.phiFY1[sc, t, s]) * mn.phiFY1[sc, t, s] + 
                              (1 - step(cap - mn.phiFY1[sc, t, s])) * cap
        mn.phiA[sc, t, s] <- step(cap - mn.phiA1[sc, t, s]) * mn.phiA1[sc, t, s] + 
                              (1 - step(cap - mn.phiA1[sc, t, s])) * cap
        mn.phiB[sc, t, s] <- step(cap - mn.phiB1[sc, t, s]) * mn.phiB1[sc, t, s] + 
                              (1 - step(cap - mn.phiB1[sc, t, s])) * cap   
  
        #Recruitment
        logit( mn.psiFYB[sc, t, s] ) <- eta[4, s, t] +  
          lmus[4, s]   # first year to breeder
        logit( mn.psiAB[sc, t, s] ) <- eta[5, s, t] +  
          lmus[5, s] # nonbreeder to breeder
        logit( mn.psiBA[sc, t, s] ) <- 
          lmus[6, s]  # breeder to nonbreeder
        #Re-encounter
        logit( mn.pA[sc, t, s] ) <- eta[7, s, t] + 
          lmus[7, s]  # resight of nonbreeders
        logit( mn.pB[sc, t, s] ) <- eta[8, s, t] +  
          lmus[8, s]  # resight of breeders
        } # s
      } # t
    } # sc
    
    # Define state-transition and observation matrices
    for (i in 1:nind){  
      # Define probabilities of state S(t+1) given S(t)
      for (t in first[i]:(nyr-1)){
        ps[1,i,t,1] <- 0
        ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t]) 
        ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
        ps[1,i,t,4] <- (1-phiFY[i,t])
        
        ps[2,i,t,1] <- 0
        ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
        ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
        ps[2,i,t,4] <- (1-phiA[i,t])
        
        ps[3,i,t,1] <- 0
        ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
        ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
        ps[3,i,t,4] <- (1-phiB[i,t])
        
        ps[4,i,t,1] <- 0
        ps[4,i,t,2] <- 0
        ps[4,i,t,3] <- 0
        ps[4,i,t,4] <- 1
        
        # Define probabilities of O(t) given S(t)
        po[1,i,t,1] <- 1 
        po[1,i,t,2] <- 0
        po[1,i,t,3] <- 0
        po[1,i,t,4] <- 0
        
        po[2,i,t,1] <- 0
        po[2,i,t,2] <- pA[i,t]
        po[2,i,t,3] <- 0
        po[2,i,t,4] <- (1-pA[i,t])
        
        po[3,i,t,1] <- 0
        po[3,i,t,2] <- 0
        po[3,i,t,3] <- pB[i,t]
        po[3,i,t,4] <- (1-pB[i,t])
        
        po[4,i,t,1] <- 0
        po[4,i,t,2] <- 0
        po[4,i,t,3] <- 0
        po[4,i,t,4] <- 1
      } #t
    } #i
    
    # Likelihood 
    for (i in 1:nind){
      # Define latent state at first capture
      z[i,first[i]] <- y[i,first[i]]
      for (t in (first[i]+1):nyr){
        # State process: draw S(t) given S(t-1)
        z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
        # Observation process: draw O(t) given S(t)
        y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
      } #t
    } #i
  } )

#**********************
#* Function to run model in NIMBLE
#**********************
run_pva <- function(info, datl, constl, code){
  library('nimble')
  library('coda')
  # function for multivariate normal
  uppertri_mult_diag <- nimbleFunction(
    run = function(mat = double(2), vec = double(1)) {
      returnType(double(2))
      p <- length(vec)
      out <- matrix(nrow = p, ncol = p, init = FALSE)
      for(i in 1:p){
        out[ ,i] <- mat[ ,i] * vec[i]
      }
      return(out)
    })
  assign('uppertri_mult_diag', uppertri_mult_diag, envir = .GlobalEnv)
  
  params <- c(
    # pop growth 
    "lambda",
    # prod
    "lmu.prod", "gamma", "rr", "mn.prod", 
    # survival 
    "mus", "lmus", "betas", "deltas",
    # abundance
    "NB", "NF", "NFY", "N", "NAD", "Ntot",
    "r",
    # error terms
    "eta", "sds2", "Ustar2", "U2", "R2",
    # yearly summaries
    'mn.phiFY', 'mn.phiA', 'mn.phiB',
    'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
    'mn.pA', 'mn.pB',
    # pva
    "extinct", "extinctAD", "extinctB"
  )
  #n.chains=1; n.thin=1; n.iter=50; n.burnin=25
  n.chains=1; n.thin=400; n.iter=1000000; n.burnin=600000
  
  mod <- nimbleModel(code, 
                     constants = constl, 
                     data = datl, 
                     inits = info$inits, 
                     buildDerivs = FALSE,
                     calculate=T 
  ) 
  
  cmod <- compileNimble(mod, showCompilerOutput = TRUE)
  confhmc <- configureMCMC(mod)
  confhmc$setMonitors(params)
  hmc <- buildMCMC(confhmc)
  chmc <- compileNimble(hmc, project = mod, 
                        resetFunctions = TRUE,
                        showCompilerOutput = TRUE)
  
  post <- runMCMC(chmc,
                  niter = n.iter, 
                  nburnin = n.burnin,
                  nchains = n.chains,
                  thin = n.thin,
                  samplesAsCodaMCMC = T, 
                  setSeed = info$seed)
  
  return(post)
} # run_pva function end

#*****************
#* Run chains in parallel
#*****************
this_cluster <- makeCluster(cpus)
post <- parLapply(cl = this_cluster, 
                  X = par_info_pva[1:cpus], 
                  fun = run_pva, 
                  datl = datl, 
                  constl = constl, 
                  code = mycode)
stopCluster(this_cluster)
save(post, mycode, seeds, cpus,
     file="/bsuscratch/brianrolek/riha_ipm/outputs/pva_longrun.rdata")
# save(post, mycode, seeds, cpus,
#      file="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//outputs//pva_survival_longrun.rdata")

2. Figures and tables

2.1 Plot model estimates

library ('MCMCvis')
library ('coda')
library ('ggplot2')
library('reshape2')
library ("tidybayes")
library ('bayestestR')
library ("ggpubr")
library("viridis")
library ("HDInterval")
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_longrun.rdata")
load("data/data.rdata")
out <- lapply(post, as.mcmc) 
outp <- MCMCpstr(out, type="chains")

#********************
#* Calculate summary stats
#********************
#* Abundance at both sites
Nall <- apply(outp$Ntot, c(1,3) , sum)
l.all.post <- array(NA, dim=c(constl$nyr-1, ncol(Nall)) ) 
for (t in 2:13){
l.all.post[t-1, ]  <- Nall[t,] / Nall[t-1,]
}
l.mn <- apply(l.all.post, 2, mean)
# Average annual pop growth
median(l.mn)
## [1] 1.020389
HDInterval::hdi(l.mn)
##     lower     upper 
## 0.9811828 1.0491875 
## attr(,"credMass")
## [1] 0.95
# 2011 to 2023 pop growth, growth between years
l.all.post2  <- Nall[13,] / Nall[1,]
median(l.all.post2)
## [1] 1.246667
HDInterval::hdi(l.all.post2)
##     lower     upper 
## 0.7111111 1.6995885 
## attr(,"credMass")
## [1] 0.95
# Population growth rates
mn.lambda <- apply(outp$lambda, c(2,3), mean)
(md.lambda <- apply(mn.lambda, 1, median))
## [1] 1.006388 1.298667
(hdi95.lambda <- apply(mn.lambda, 1, HDInterval::hdi))
##            [,1]     [,2]
## lower 0.9669408 1.206648
## upper 1.0360502 1.420549
# Overall averages
mn.mus <- apply(outp$mus, c(1,3), mean)
(md.mus <- apply(mn.mus, 1, median)) |> round(2)
## [1] 0.31 0.94 0.89 0.12 0.43 0.07 0.21 0.77
(hdi.mus <- apply(mn.mus, 1, HDInterval::hdi)) |> round(2)
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.17 0.82 0.82 0.02 0.28 0.03 0.04 0.53
## upper 0.47 1.00 0.96 0.31 0.55 0.12 0.43 0.96
# Compare between sites
df.sites <- data.frame(mus=1:8, pd=rep(NA,8), greater=NA)
for (i in 1:8){
  first <- median(outp$mus[i, 1, ]) > median(outp$mus[i, 2, ])
  df.sites$greater[i] <- ifelse(first, "LH", "PC")
  if (first){
    diff <- outp$mus[i, 1, ]-outp$mus[i, 2, ]
    df.sites$pd[i] <- mean(diff>0) |> round(2) 
  } else{
    diff <- outp$mus[i, 2, ]-outp$mus[i, 1, ]
    df.sites$pd[i] <- mean(diff>0) |> round(2)
  }
}
df.sites
(md.mus <- apply(outp$mus, c(1,2), median)) |> round(2)
##      [,1] [,2]
## [1,] 0.35 0.26
## [2,] 0.97 0.92
## [3,] 0.91 0.87
## [4,] 0.14 0.07
## [5,] 0.12 0.73
## [6,] 0.12 0.03
## [7,] 0.06 0.33
## [8,] 0.68 0.88
(hdi.mus <- apply(outp$mus, c(1,2), HDInterval::hdi)) |> round(2)
## , , 1
## 
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.17 0.89 0.84 0.01 0.05 0.05 0.00 0.27
## upper 0.57 1.00 0.98 0.38 0.19 0.20 0.25 0.97
## 
## , , 2
## 
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.08 0.68 0.74 0.00 0.45 0.00 0.04 0.57
## upper 0.48 1.00 0.97 0.38 0.95 0.08 0.73 1.00
# productivity
(outp$lmu.prod[1,]-outp$lmu.prod[2,]) |> pd()
## Probability of Direction: 1.00
# Compare survival, recruitment, and detection rates
# among life stages within each site
df.comp <- data.frame(param1 = c('phiA', 'phiA', 'phiB', 'psiA', 'pB'), 
                 param2 = c('phiB', 'phiFY', 'phiFY', 'psiFY', 'pA'), 
                 pd_LH = NA, pd_PC = NA)
## Los Haitises, survival
df.comp$pd_LH[1] <- (outp$lmus[2,1,]-outp$lmus[3,1,]) |> pd() |> round(2)
df.comp$pd_LH[2] <- (outp$lmus[2,1,]-outp$lmus[1,1,]) |> pd() |> round(2)
df.comp$pd_LH[3] <- (outp$lmus[3,1,]-outp$lmus[1,1,]) |> pd() |> round(2)
## Los Haitises, recruitment rates
df.comp$pd_LH[4] <- (outp$lmus[5,1,]-outp$lmus[4,1,]) |> pd() |> round(2)
## Los Haitises, resight probability
df.comp$pd_LH[5] <- (outp$lmus[8,1,]-outp$lmus[7,1,]) |> pd() |> round(2)
## Punta Canas, survival
df.comp$pd_PC[1] <- (outp$lmus[2,2,]-outp$lmus[3,2,]) |> pd() |> round(2)
df.comp$pd_PC[2] <- (outp$lmus[2,2,]-outp$lmus[1,2,]) |> pd() |> round(2)
df.comp$pd_PC[3] <- (outp$lmus[3,2,]-outp$lmus[1,2,]) |> pd() |> round(2)
## Punta Cana, recruitment rates
df.comp$pd_PC[4] <- (outp$lmus[5,2,]-outp$lmus[4,2,]) |> pd() |> round(2)
## Punta Cana, resight probability
df.comp$pd_PC[5] <- (outp$lmus[8,2,]-outp$lmus[7,2,]) |> pd() |> round(2)
df.comp 
#*******************
#* Plot IPM results
#*******************
#Plot population size
lNall <- melt(Nall)
colnames(lNall) <- c("Time", "Iter", "Abundance")
lNall$Site.ind <- 3; lNall$Site <- "Both"
lNall <- lNall[,c(1,4,2,3,5)]

Ntot <- outp$Ntot[1:13,,] # thin by 10 to speed up plotting
lN <- melt(Ntot)
colnames(lN) <- c("Time", "Site.ind", "Iter", "Abundance")
lN$Site <- ifelse(lN$Site.ind==1, "Los Haitises", "Punta Cana")
lN <- rbind(lN, lNall)

lN$Year <- lN$Time + 2010

med_hdis <- function(x, cm){
            df <- data.frame(y=median(x),
                             ymin=HDInterval::hdi(x, credMass=cm)[1],
                             ymax=HDInterval::hdi(x, credMass=cm)[2] )
}

counts.tot <- datl$countsAdults+datl$countsFY+constl$hacked.counts[,-3]
df.counts <- melt(counts.tot)
colnames(df.counts) <- c("Year", "Site_ab", "Count")
df.counts$Site <- ifelse(df.counts$Site_ab=="LHNP", "Los Haitises", "Punta Cana")

df.counts.all <- data.frame(Year=2011:2023, Site_ab="Both", Count=counts.tot[,1] + counts.tot[,2], Site="Both" )
df.counts <- rbind(df.counts, df.counts.all)

p1 <- ggplot() + theme_minimal(base_size=14) + 
  geom_line(data=lN, aes(x=Year, y=Abundance, group=Iter), 
            color="gray30", linewidth=0.1, alpha=0.01 ) +
  stat_summary(data=lN, aes(x=Year, y=Abundance), 
               geom="errorbar", width=0, 
               fun.data =med_hdis, fun.args=list(cm=0.95) ) +
  stat_summary(data=lN, aes(x=Year, y=Abundance), 
               geom="errorbar", width=0, linewidth=1,
               fun.data =med_hdis, fun.args=list(cm=0.80) ) +
  stat_summary(data=lN, aes(x=Year, y=Abundance), 
               geom="line", linewidth=1,
               fun.data =function(x){data.frame(y=median(x))} ) +
  geom_line(data=df.counts, aes(x=Year, y=Count), 
            col="black", linewidth=1, linetype=2 ) +
  facet_wrap("Site", scales="free_y") +
  xlab("Year") + ylab("Number")
p1

# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Abundance and counts.tiff",
#        plot=p1,
#        device="tiff",
#        width=9, height=4, dpi=300)
# Plot life stage structure
mdFY <-  apply(outp$NFY[1:13,,], c(1,2), median) 
mdB <-  apply(outp$NB[1:13,,], c(1,2), median) 
mdF <-  apply(outp$NF[1:13,,], c(1,2), median) 
lFY <- melt(mdFY)
lB <- melt(mdB)
lF <- melt(mdF)
lFY$Stage <- "First-year"
lB$Stage <- "Breeder"
lF$Stage <- "Nonbreeder"
ldat <- rbind(lFY, lB, lF)
colnames(ldat)[1:3] <- c("Year", "Sitenum", "Number") 
ldat$Site <- ifelse(ldat$Sitenum==1, "Los Haitises", "Punta Cana")
ldat$year <- ldat$Year+2010

# Use median number of females in each stage
# to plot an approximate population structure
p3 <- ggplot(ldat, aes(fill=Stage, y=as.numeric(Number), x=year)) + 
  geom_bar(position="fill", stat="identity") +
  ylab("Proportion of population") + xlab("Year") +
  facet_wrap("Site") + 
  theme_minimal(base_size=14)+ theme(legend.position="top") +
  scale_fill_viridis_d(option = "mako")

p4 <- ggplot(ldat, aes(fill=Stage, y=as.numeric(Number), x=year)) + 
  geom_bar(position="stack", stat="identity") +
  ylab("Median number of females") + xlab("Year") +
  facet_wrap("Site", scales = "free_y") + 
  theme_minimal(base_size=14 ) + theme(legend.position="top") + 
  scale_fill_viridis_d(option = "mako") 

# plot with uncertainty
lFY <- melt(outp$NFY)
lB <- melt(outp$NB)
lF <- melt(outp$NF)
colnames(lFY) <- colnames(lB) <- colnames(lF) <- c("Time", "Site", "Iteration", "Abundance")
lFY$Year <- lFY$Time + 2010
lB$Year <- lB$Time + 2010
lF$Year <- lF$Time + 2010
lFY$Stage <- "First-year"
lB$Stage <- "Breeder"
lF$Stage <- "Nonbreeder"
lALL <- rbind(lFY, lB, lF)
lALL$Site <- ifelse(lALL$Site==1, "Los Haitises", "Punta Cana")

mdFY <-  apply(outp$NFY[1:13,,], c(1,2), median) 
mdB <-  apply(outp$NB[1:13,,], c(1,2), median) 
mdF <-  apply(outp$NF[1:13,,], c(1,2), median) 
hdiFY <-  apply(outp$NFY[1:13,,], c(1,2), HDInterval::hdi) 
hdiB <-  apply(outp$NB[1:13,,], c(1,2), HDInterval::hdi) 
hdiF <-  apply(outp$NF[1:13,,], c(1,2), HDInterval::hdi) 
medFY <- melt(mdFY)
medB <- melt(mdB)
medF <- melt(mdF)
hdisFY <- melt(hdiFY)
hdisB <- melt(hdiB)
hdisF <- melt(hdiF)
dfstages <- rbind(medFY, medB, medF)
dfhdis <- rbind(hdisFY, hdisB, hdisF)
dfstages$LHDI <- dfhdis[dfhdis$Var1=="lower",]
dfstages$UHDI <- dfhdis[dfhdis$Var1=="upper",]
cols <- mako(3)

med_hdis <- function(x, cm){
  df <- data.frame(y=median(x),
                   ymin=HDInterval::hdi(x, credMass=cm)[1],
                   ymax=HDInterval::hdi(x, credMass=cm)[2] )
}

p5 <- ggplot() +  theme_minimal(base_size=14) +
  theme(legend.position="none") +
  geom_line(data=lALL, aes(x=Year, y=Abundance, 
                           group=Iteration, color=Stage,
                           alpha=Stage), 
            linewidth=0.1) +
  stat_summary(data=lALL, aes(x=Year, y=Abundance),
               geom="errorbar", width=0,
               fun.data =med_hdis, fun.args=list(cm=0.95) ) +
  stat_summary(data=lALL, aes(x=Year, y=Abundance),
               geom="errorbar", width=0, linewidth=1,
               fun.data =med_hdis, fun.args=list(cm=0.80) ) +
  stat_summary(data=lALL, aes(x=Year, y=Abundance),
               geom="line", linewidth=1,
               fun.data =function(x){data.frame(y=median(x))} ) +
  scale_color_manual( values = c("First-year"=cols[2], "Breeder"=cols[1], "Nonbreeder"=cols[3])  ) +
  scale_alpha_manual( values = c("First-year"=0.015, "Breeder"=0.0075, "Nonbreeder"=0.15) ) +
  facet_wrap(Stage~Site, scales="free_y", 
             shrink=TRUE, ncol=2) +
  xlab("Year") + ylab("Number of females") +
  theme(
    strip.background = element_blank(),
    strip.text.x = element_blank()
  )

p35 <- ggarrange(p3, p5, ncol=1, nrow=2)
p35

# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Stage structure_abundance.tiff",
#        plot=p35,
#        device="tiff",
#        width=6, height=8, dpi=300)
# plot survival, recruitment, and detection
mus.mat <- outp$mus
dimnames(mus.mat)[[1]] <- c("FY", 
                            "A\nSurvival", "B", 
                            "FY:B",
                            "A:B\nRecruitment", 
                            "B:A",
                            "A", "B\nDetection")
dimnames(mus.mat)[[2]] <- c("Los Haitises", "Punta Cana")
lmus <- melt(mus.mat)
colnames(lmus) <- c("Parameter", "Site", "Iter", "value")

# p <- list()
# p[[1]] <- apply(outp$mn.phiFY, c(2, 3), mean)
# p[[2]] <- apply(outp$mn.phiA, c(2, 3), mean)
# p[[3]] <- apply(outp$mn.phiB, c(2, 3), mean)
# p[[4]] <- apply(outp$mn.psiFYB, c(2, 3), mean)
# p[[5]] <- apply(outp$mn.psiAB, c(2, 3), mean)
# p[[6]] <- apply(outp$mn.psiBA, c(2, 3), mean)
# p[[7]] <- apply(outp$mn.pA, c(2, 3), mean)
# p[[8]] <- apply(outp$mn.pB, c(2, 3), mean)
# pl <- lapply(p, melt)
# for (i in 1:8){
#   nms <- c('FY', 'A', 'B', 'FY:B', 'A:B', 'B:A', 'pA', 'pB') 
#   pl[[i]]$nm <- nms[i]
# } 
# pl <- do.call(rbind, pl)
# colnames(pl) <- c('Site', 'Iter', 'value', 'nm')
# pl$nm <- factor(pl$nm, levels=c('FY', 'A', 'B', 'FY:B', 'A:B', 'B:A', 'pA', 'pB'))

p6 <- ggplot(data= lmus, aes(x=value, y=Parameter )) +
  stat_halfeye(point_interval=median_hdi, .width = c(.95, .80)) +
  coord_flip() +
  facet_wrap("Site") +
  xlab("Probability") + ylab("Parameter") +
  theme_bw(base_size=14) +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()) 

p6

# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Survival Recruitment Detection.tiff",
#        plot=p5,
#        device="tiff",
#        width=8, height=4, dpi=300)

# print estimates
mds <- tapply(lmus$value, list(lmus$Parameter, lmus$Site), median) |> round(2)
hdi_LH <- tapply(lmus$value, list(lmus$Parameter, lmus$Site), HDInterval::hdi)[,1]
hdi_LH2 <- do.call(rbind, hdi_LH) |> round(2)                                                                  
hdi_PC <- tapply(lmus$value, list(lmus$Parameter, lmus$Site), HDInterval::hdi)[,2]
hdi_PC2 <- do.call(rbind, hdi_PC) |> round(2) 
df <- data.frame(param= c("phiFY", "phiA", "phiB", "psiFY:B", "psiA:B", "psiB:A", "pA", "pB"),
  md_LH=mds[,1], lhdi_LH=hdi_LH2[,1], uhdi_LH=hdi_LH2[,2], 
  md_PC=mds[,2], lhdi_PC=hdi_PC2[,1], uhdi_PC=hdi_PC2[,2])
df
# Plot predicted survival, recruitment, and detection
# with effects from translocation and hacking
# Birds were only hacked from LHNP to PC here
# so we only predict values for PC

betas.pd <- apply(outp$betas, 1, pd)
ind <- which (betas.pd>0.975) # which are significant
pred.mus <- array(NA, dim=c(dim(outp$mus)))
for (m in ind){
  for (tr in 1:2){
    pred.mus[m,tr,] <- plogis( outp$lmus[m,2,] + outp$betas[m,]*c(0,1)[tr] ) 
  }}
# treatment, mu, site, iter
mus.md <- apply(pred.mus, c(1,2), median)
mus.HDI80 <- apply(pred.mus, c(1,2), HDInterval::hdi, credMass=0.8)
mus.HDI95 <- apply(pred.mus, c(1,2), HDInterval::hdi)
# tiff(height=4, width=6, units="in", res=300,
#      filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Translocation effects.tiff")
par(mar=c(3,4,4,1), mfrow=c(1,1))
for (tr in 1:2){
  if (tr==1){
    plot(1:4, c(0,0,1,1), 
         type="n",
         pch=c(16, 17)[tr], cex=2,
         ylim=c(0,1),
         xlim=c(0.5, 4.5), 
         cex.axis=1, cex.lab=1.25,
         ylab="Probability", xlab="",
         main="",
         xaxt="n", yaxt="n")
  abline(h= seq(0,1,by=0.1), col="gray90")
  abline(v= seq(0.5,6,by=1), col="gray90")
  points(1:4+c(-0.1,0.1)[tr], mus.md[ind,tr], 
         pch=c(16, 17)[tr], cex=1.5)
  }else{
    points(1:4+c(-0.1,0.1)[tr], mus.md[ind,tr], 
           pch=c(16, 15)[tr], cex=1.5)
  }}

axis(1, at=1:4, cex.axis=0.85,
     labels=c("First-year\nSurvival", 
              "Breeder\nSurvival", 
              "Nonbreeder\n to Breeder\nRecruitment", 
              "Nonbreeder\nDetection"), las=1, padj=0.5)
axis(2, at=seq(0,1, by=0.1), cex.axis=1,
     labels=c(0, NA,NA,NA,NA,0.5, NA,NA,NA,NA, 1))
for (m in 1:length(ind)){
  for (tr in 1:2){
    lines(rep(c(1:4)[m] + c(-0.1,0.1)[tr],2), mus.HDI95[1:2,ind[m],tr], lwd=2)
    lines(rep(c(1:4)[m] + c(-0.1,0.1)[tr],2), mus.HDI80[1:2,ind[m],tr], lwd=4)
  }}
legend(x=2.5,y=1.25, pch=c(16,15), pt.cex=1.5, cex=1, xpd=NA, horiz=T, 
       legend=c("Not translocated", "Translocated" ),
       xjust=0.5)

#dev.off()
# Plot fecundity in response to nest treatments
f <- exp(outp$lmu.prod)
f.md <- apply(f, 1, median)
f.HDI80 <- apply(f, 1, HDInterval::hdi, credMass=0.8)
f.HDI95 <- apply(f, 1, HDInterval::hdi)

# Calculate treatment effects on fecundity
f.pred <- array(NA, dim=dim(outp$lmu.prod))
for (s in 1:2){
  f.pred[s,] <- exp(outp$lmu.prod[s,] + outp$gamma[,1])
} # s
f2.md <- apply(f.pred, 1, median)
f2.HDI80 <- apply(f.pred, 1, HDInterval::hdi, credMass=0.8)
f2.HDI95 <- apply(f.pred, 1, HDInterval::hdi)

# tiff(height=4, width=6, units="in", res=300,
#      filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Fecundity_treatment.tiff")
par(mar=c(3,5,3,1), mfrow=c(1,1))
plot(c(1, 2)-0.1, f.md, 
     ylim=c(0, 1.55),
     xlim=c(0.5, 2.5), 
     ylab="Productivity", xlab="",
     main="",
     cex.axis=1.5, cex.lab=1.5,
     xaxt="n", yaxt="n", type="n")
points(c(1, 2)-0.1, f.md, pch=16, cex=1.5)
abline(h=seq(0, 1.6, by=0.1), col="gray90")
abline(v=1.5, col="black")
axis(1, at=c(1,2), cex.axis=1,
     labels=c("Los Haitises\nNational Park","Punta Cana"),
     padj=0.5)
axis(2, at=seq(0, 1.6, by=0.1), 
     cex.axis=1,
     labels=c("0", NA, NA, NA, NA, "0.5", NA, NA, NA, NA, "1.0", NA, NA, NA, NA, "1.5", NA))
lines(c(1,1)-0.1, f.HDI95[,1], lwd=2)
lines(c(2,2)-0.1, f.HDI95[,2], lwd=2)
lines(c(1,1)-0.1, f.HDI80[,1], lwd=4)
lines(c(2,2)-0.1, f.HDI80[,2], lwd=4)

points(c(1.1, 2.1), f2.md, 
       pch=18, cex=2.5)

lines(c(1,1) +0.1, f2.HDI95[,1], lwd=3)
lines(c(2,2) +0.1, f2.HDI95[,2], lwd=3)
lines(c(1,1) +0.1, f2.HDI80[,1], lwd=6)
lines(c(2,2) +0.1, f2.HDI80[,2], lwd=6)
legend(x=1.5,y=1.9, pch=c(16,18), pt.cex=c(1.5,2.5), cex=1,
       legend=c("Untreated", "Treated" ), 
       horiz=TRUE, xjust=0.5, xpd=NA)

#dev.off()

# Calculate probability of direction
f.diff <- f[1,]-f[2,]
pd(f.diff)
## Probability of Direction: 1.00
f.md 
## lmu.prod[1] lmu.prod[2] 
##   0.3743717   0.2265712
f.HDI95 
##       lmu.prod[1] lmu.prod[2]
## lower   0.2952690   0.1640874
## upper   0.4527469   0.2986755
f.diff2 <- f.pred[1,]-f.pred[2,]
pd(f.diff2)
## Probability of Direction: 1.00
f2.md 
## [1] 1.2559463 0.7601034
f2.HDI95
##            [,1]     [,2]
## lower 0.9905714 0.550482
## upper 1.5188800 1.002000
median(f.pred[1,]/f[1,])
## [1] 3.354811
median(f.pred[2,]/f[2,])
## [1] 3.354811
#*******************
# Correlations between demographic rates
#*******************
pdR2 <- apply(outp$R2, c(1,2), pd) |> round(2)
max(pdR2[pdR2<1])
## [1] 0.82
apply(outp$R2, c(1,2), median) |> round(2)
##        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9]
##  [1,]  1.00 -0.02  0.14  0.11 -0.07  0.00  0.08  0.07  0.05
##  [2,] -0.02  1.00 -0.02 -0.03  0.02  0.00 -0.07 -0.02  0.00
##  [3,]  0.14 -0.02  1.00 -0.05 -0.11  0.01  0.20 -0.06 -0.04
##  [4,]  0.11 -0.03 -0.05  1.00 -0.11  0.00  0.23 -0.22 -0.17
##  [5,] -0.07  0.02 -0.11 -0.11  1.00 -0.01 -0.08 -0.02  0.08
##  [6,]  0.00  0.00  0.01  0.00 -0.01  1.00  0.01  0.00  0.00
##  [7,]  0.08 -0.07  0.20  0.23 -0.08  0.01  1.00 -0.22 -0.21
##  [8,]  0.07 -0.02 -0.06 -0.22 -0.02  0.00 -0.22  1.00  0.04
##  [9,]  0.05  0.00 -0.04 -0.17  0.08  0.00 -0.21  0.04  1.00
apply(outp$R2, c(1,2), HDInterval::hdi)[1,,] |> round(2)
##        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9]
##  [1,]  1.00 -0.64 -0.45 -0.50 -0.64 -0.55 -0.46 -0.45 -0.52
##  [2,] -0.64  1.00 -0.61 -0.62 -0.57 -0.57 -0.63 -0.59 -0.59
##  [3,] -0.45 -0.61  1.00 -0.63 -0.66 -0.62 -0.35 -0.57 -0.54
##  [4,] -0.50 -0.62 -0.63  1.00 -0.68 -0.57 -0.34 -0.71 -0.71
##  [5,] -0.64 -0.57 -0.66 -0.68  1.00 -0.62 -0.61 -0.57 -0.49
##  [6,] -0.55 -0.57 -0.62 -0.57 -0.62  1.00 -0.60 -0.62 -0.63
##  [7,] -0.46 -0.63 -0.35 -0.34 -0.61 -0.60  1.00 -0.63 -0.65
##  [8,] -0.45 -0.59 -0.57 -0.71 -0.57 -0.62 -0.63  1.00 -0.45
##  [9,] -0.52 -0.59 -0.54 -0.71 -0.49 -0.63 -0.65 -0.45  1.00
apply(outp$R2, c(1,2), HDInterval::hdi)[2,,] |> round(2)
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
##  [1,] 1.00 0.56 0.67 0.62 0.51 0.62 0.62 0.59 0.60
##  [2,] 0.56 1.00 0.55 0.58 0.59 0.60 0.53 0.58 0.59
##  [3,] 0.67 0.55 1.00 0.53 0.48 0.58 0.69 0.45 0.44
##  [4,] 0.62 0.58 0.53 1.00 0.44 0.60 0.71 0.31 0.40
##  [5,] 0.51 0.59 0.48 0.44 1.00 0.56 0.48 0.52 0.61
##  [6,] 0.62 0.60 0.58 0.60 0.56 1.00 0.57 0.59 0.59
##  [7,] 0.62 0.53 0.69 0.71 0.48 0.57 1.00 0.23 0.33
##  [8,] 0.59 0.58 0.45 0.31 0.52 0.59 0.23 1.00 0.50
##  [9,] 0.60 0.59 0.44 0.40 0.61 0.59 0.33 0.50 1.00
#*******************
# Correlations between demographics and population growth rates
#*******************
lam.m <- apply(outp$lambda[1:12,,], c(1,2), median)
lam.hdi <- apply(outp$lambda[1:12,,], c(1,2), HDInterval::hdi)
par(mfrow=c(1,2))
plot(2012:2023, lam.m[,1], type="b", pch=1, 
     ylab="Population growth rate", xlab="Year", 
     ylim=c(min(lam.hdi[,,1]), max(lam.hdi[,,1])),
     main="Los Haitises")
abline(h=1, lty=2)
segments(x0=2012:2023, x1=2012:2023, 
         y0 = lam.hdi[1,,1], y1= lam.hdi[2,,1])

plot(2012:2023, lam.m[,2], type="b", pch=1, lty=1,
     ylab="Population growth rate", xlab="Year",
     ylim=c(min(lam.hdi[,,2]), max(lam.hdi[,,2])),
     main="Punta Cana")
abline(h=1, lty=2)
segments(x0=2012:2023, x1=2012:2023, 
         y0 = lam.hdi[1,,2], y1= lam.hdi[2,,2])

# create a function to plot correlations
# between demographics and population growth rates
plot.cor <- function (lam.post, x.post, x.lab, ind.x=1:12, site=1, yaxt="b"){
  # calculate the correlation coefficient 
  # over each iteration to propagate error
  lambda <- lam.post[1:12, site,]
  x <- x.post[ind.x, site, ]
  df <- data.frame(ats1=c(0.8, 0.9, 1.0, 1.1, 1.2, 1.3),
                   ats2=c(1.0, 1.5, 2.0, 2.5, 3.0, 3.5),
                   labs1=c("0.8", NA, "1.0", NA, "1.2", NA),
                   labs2=c("1.0", NA, "2.0", NA, "3.0", NA)
  )
  if(site==1){ df <- df[, c(1,3)]} else{
    df <- df[, c(2,4)]
  }
  lwd <- 1
  
  cor.post <- array(NA, dim=dim(lambda)[-1])
    for (i in 1:dim(lambda)[[2]]){
      cor.post[i] <- cor(lambda[,i], x[,i])
    }
  cor.df <- data.frame(median= median(cor.post) |> round(digits=2),
                       ldi=HDInterval::hdi(cor.post)[1] |> round(digits=2),
                       hdi=HDInterval::hdi(cor.post)[2] |> round(digits=2),
                       pd = pd(cor.post) |> round(digits=2))
  
  lam.m <- apply(lambda, 1, median)
  lam.hdi <- apply(lambda, 1, HDInterval::hdi)
  x.m <- apply(x, 1, median)
  x.hdi <- apply(x, 1, HDInterval::hdi)
    x.lims <- c(min(x.hdi[,]), max(x.hdi[,]))
    y.lims <- c(min(lam.hdi[,]), max(lam.hdi[,]))
    if(yaxt=="n"){
    plot(x.m, lam.m[1:12], 
         xlim= x.lims,
         ylim= y.lims,
         type="n", 
         ylab="", xlab=x.lab, cex.lab=1.5, 
         main=c("", "")[site],
         yaxt=yaxt, xaxt="n")
      axis(2, at=df[,1], labels=c(rep(NA, 5)))
    } else {
      plot(x.m, lam.m[1:12], 
           xlim= x.lims,
           ylim= y.lims,
           type="n", 
           ylab="", xlab=x.lab, cex.lab=1.5,  
           main=c("", "")[site], 
           yaxt="n", xaxt="n")
      axis(2, at=df[,1], labels=df[,2], las=1, cex.axis=1.5)
    }
    axis(1, cex.axis=1.5)
    points(x.m, lam.m, pch=1)
    segments(x0=x.hdi[1,], x1=x.hdi[2,], 
             y0 = lam.m, y1= lam.m, lwd=lwd)
    segments(x0=x.m, x1=x.m, 
             y0 = lam.hdi[1,], y1= lam.hdi[2,], lwd=lwd)
    xs <- c(x.lims[1], x.lims[1]+(x.lims[2]-x.lims[1]*1.5))
    ys <- c( (y.lims[2]-y.lims[1])+y.lims[1], (y.lims[2]-y.lims[1])+y.lims[1], 
             (y.lims[2]-y.lims[1])*0.6+y.lims[1], (y.lims[2]-y.lims[1])*0.6+y.lims[1])
    polygon(x=c(xs, rev(xs)), y=ys, col=alpha("white", 0.7), border=NA )
    text(x = x.lims[1], y = (y.lims[2]-y.lims[1])*0.9+y.lims[1], 
         paste("r = ", cor.df$median, " (",cor.df$ldi,", ", cor.df$hdi, ")", sep=""), 
         pos = 4, font = 3, cex = 1.5)
    text(x = x.lims[1], y = (y.lims[2]-y.lims[1])*0.7+y.lims[1], paste("P(r>0) = ", cor.df$pd, sep=""), 
         pos = 4, font = 3, cex = 1.5)
    
  }

# Plor correlations between population grrowth rates
# and demographics. 
# "r" is a correlation coefficient and represents
# the magnitude of the correlation
# P(r>0) is the probability of direction (similar to p-values)
# that is, the probability that an effect exists
# tiff(height=8, width=8, units="in", res=300,
#      filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//PopGrowthand Demographics1.tiff")
# par(mfrow=c(4,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
par(mfrow=c(2,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
plot.cor(outp$lambda, outp$mn.prod, x.lab="Fecundity", ind.x=2:13)
plot.cor(outp$lambda, outp$mn.phiFY, x.lab="First-year Survival", yaxt="n")
plot.cor(outp$lambda, outp$mn.phiA, x.lab="Nonbreeder Survival", yaxt="n")
plot.cor(outp$lambda, outp$mn.phiB, x.lab="Breeder Survival")
plot.cor(outp$lambda, outp$mn.psiFYB, x.lab="First-year to Breeder", yaxt="n")
plot.cor(outp$lambda, outp$mn.psiAB, x.lab="Nonbreeder to Breeder", yaxt="n")
mtext("Population growth rate", side=2, outer=TRUE, line=2.5, cex=2)
mtext("Los Haitises", side=3, outer=TRUE, cex=2)

#dev.off()

# tiff(height=4, width=8, units="in", res=300,
#      filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//PopGrowthand Demographics2.tiff")
par(mfrow=c(2,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
plot.cor(outp$lambda, outp$mn.prod, x.lab="Fecundity", ind.x=2:13, site=2)
plot.cor(outp$lambda, outp$mn.phiFY, x.lab="First-year Survival", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.phiA, x.lab="Nonbreeder Survival", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.phiB, x.lab="Breeder Survival", site=2)
plot.cor(outp$lambda, outp$mn.psiFYB, x.lab="First-year to Breeder", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.psiAB, x.lab="Nonbreeder to Breeder", yaxt="n", site=2)
mtext("Population growth rate", side=2, outer=TRUE, line=2.5, cex=2)
mtext("Punta Cana", side=3, outer=TRUE, cex=2)

#dev.off()
#*******************
#* Plot PVA results
#*******************
library ('viridis')
library ("tidybayes")
library ('coda')
library ('MCMCvis')
library ('reshape2')
library('ggplot2')
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\pva_longrun.rdata")
out <- lapply(post, as.mcmc) 
outp <- MCMCpstr(out, type="chains")
scen <- data.frame('Scenarios' = 1:36,
                   'Number translocated' = rep(c(0,5,10), each=12), 
                   'Number of nests treated' = rep( rep(c(0, 15, 30, 45), each=3), 3),
                   'Mortality reduction' = rep(c('100% Mortality', '80% Mortality', '60% Mortality'), 12) )

#*********************
#* Line plot of extinction probability over time
#*********************
pe <- apply(outp$extinct, c(1,2,3), mean)
lpe <- melt(pe)
colnames(lpe) <- c('Scenarios', 'Time', 'Site2', 'Extinct')
pe.df <- merge(lpe, scen, by='Scenarios', all.x=TRUE)
pe.df$Year <- pe.df$Time+2023
pe.df$Site <- ifelse(pe.df$Site2==1, "Los Haitises", "Punta Cana")
pe.df[(pe.df$Site2==2 & pe.df$Number.translocated>0),]$Extinct <- NA
pe.df$Mortality.reduction <- factor(pe.df$Mortality.reduction, 
                                    levels= c('100% Mortality', '80% Mortality', '60% Mortality'))

p6 <- ggplot(pe.df, aes(x=Year, y=Extinct, group=Scenarios, 
                  color = factor(Number.translocated), 
                  linetype = factor(Number.of.nests.treated))) +
  geom_line( linewidth = 1) +
  facet_grid(rows=  vars(Site), cols= vars(Mortality.reduction), 
             scales='free_y') +
  #facet_wrap(~  Site + Mortality.reduction) +
  scale_color_viridis_d(option="viridis", direction=1,
                        name="Number translocated") +
  scale_linetype(name="Number of nests treated") +
  scale_x_continuous( name="Future year",
                      breaks=c(2020, 2030, 2040, 2050, 2060, 2070),
                     labels=c('', '2030', '', '2050', '', '2070'), 
                     minor_breaks=NULL) +
  theme_minimal() + 
  ylab("Extinction probability")
p6

# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Extinction_lines.tiff",
#        plot=p6,
#        device="tiff",
#        width=6, height=4, dpi=300)

#********************
#* PVA heatmap of extinction probability projected 50 years
#********************
pe.df2 <- pe.df[pe.df$Time==50, ]

p7 <- ggplot(pe.df2, aes(x = as.factor(Number.translocated), 
                  y = as.factor(Number.of.nests.treated))) +
  geom_tile( aes(fill = Extinct ) ) + 
  geom_text(aes(label=round(Extinct,2)), color="gray70") +
  scale_fill_viridis_c(option="plasma", direction=1,
                       name = "Extinction\nprobability after\n50 years") +
  facet_wrap(~  Site + factor(Mortality.reduction, 
                              levels=c("100% Mortality", "80% Mortality", "60% Mortality"))) +
  scale_y_discrete(breaks=c(0,15,30,45), 
                   name="Number of nests treated") +
  scale_x_discrete(name="Number translocated") +
  theme_minimal( ) + theme(strip.background = element_rect(size = 0.5))
p7  

# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Extinction_50yrHeatmap.tiff",
#        plot=p7,
#        device="tiff",
#        width=6, height=4, dpi=300)
pe.df3 <- pe.df2[!is.na(pe.df2$Extinct),]
pe.df3 <- pe.df3[rev(order(pe.df3$Site2, pe.df3$Extinct)),]
pe.df3
#***************************
# Plot total Abundance
lnt1 <-  outp$Ntot[,1:13, , ] |> melt()
colnames(lnt1) <- c("Scenarios", "Time", "Site2", "Iter", "Abund")
lnt1$Year <- lnt1$Time + 2010
lnt1$Site <- factor(ifelse(lnt1$Site2==1, 'Los Haitises', 'Punta Cana'))
nt.df <- merge(lnt1, scen, by='Scenarios', all.x=TRUE)
nt.df$Mortality.reduction <- factor(nt.df$Mortality.reduction,
                                    levels=c('100% Mortality','80% Mortality', '60% Mortality' ))
nt.df[(nt.df$Site2==2 & nt.df$Number.translocated>0),]$Abund <- NA
nt.df <- nt.df[nt.df$Scenarios %in% c(1,2,3), ]

# calculate future means
lnt2 <- apply(outp$Ntot[,13:63, , ], c(1, 2, 3), median) |> melt()
colnames(lnt2) <- c("Scenarios", "Time", "Site2", "Abund")
lnt2$Year <- lnt2$Time + 2022
lnt2$Site <- factor(ifelse(lnt2$Site2==1, 'Los Haitises', 'Punta Cana'))
nt.df2 <- merge(lnt2, scen, by='Scenarios', all.x=TRUE)
nt.df2$Mortality.reduction <- factor(nt.df2$Mortality.reduction,
                                    levels=c('100% Mortality','80% Mortality', '60% Mortality' ))
nt.df2[(nt.df2$Site2==2 & nt.df2$Number.translocated>0),]$Abund <- NA

p8 <- ggplot() +
      stat_lineribbon(data= nt.df,  
                      aes(x=Year, y=Abund, group=Scenarios),
                      fill='gray40',
                      linewidth=0.5,
                      point_interval = 'median_hdci', .width = 0.95,
                       na.rm=TRUE) +
      geom_line(data= nt.df2,
                aes(x=Year, y=Abund, group=Scenarios,
                    color = factor(Number.translocated),
                    linetype = factor(Number.of.nests.treated) ),
                linewidth=0.5) +
      facet_grid(rows = vars(Site), cols = vars(Mortality.reduction),
                 scales ='free_y') +
      ylab('Number of females') +
      scale_color_viridis_d(option="viridis", direction=1,
                            name="Number translocated") +
      scale_linetype(name="Number of nests treated") +
      theme_minimal()
  
p8

# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Abundance_projections.tiff",
#        plot=p9,
#        device="tiff",
#        width=8, height=4, dpi=300)

# Time to extinction
# conditional on extinction threshold (D=X)
D <- 1
abund <- outp$Ntot[,14:(13+50),,]
dimnames(abund) <- list('Scenarios'=1:36, 
                        'Time'= 2024:(2024+49), 
                        'Site'=c("Los Haitises", "Punta Cana"), 
                        'Iter'=1:10000)
qextinct <- abund == D
ext <- qextinct[,50,,] <= D
lte <- melt( qextinct )
colnames(lte)[5] <- "Extinct"
le <- melt( ext )
colnames(le)[4] <- "EverExtinct"
lte2 <- merge(lte, le, by=c("Scenarios", "Site", "Iter")  )
lte3<- lte2[lte2$EverExtinct==TRUE, ]
lte4 <- lte3[lte3$Extinct==TRUE, ]
lte5 <- tapply(lte4$Time-2023, list(lte4$Scenarios, lte4$Site, lte4$Iter), min)
lte6 <- melt(lte5)
lte7 <- lte6[!is.na(lte6$value), ]
colnames(lte7) <- c('Scenarios', 'Site', 'Iter', 'TTE')
te <- table(lte7$TTE, lte7$Scenarios, lte7$Site) |> melt()
colnames(te) <- c('TE', 'Scen', 'Site', 'Count')
te$Scenarios <- factor(paste0("Scenario ", te$Scen), 
                       levels=paste0('Scenario ', 1:36)) 
teLH <- te[te$Site=="Los Haitises",]
tePC <- te[te$Site=="Punta Cana",]

p9 <- ggplot() +
      geom_col(data = teLH, aes(x= TE, y=Count)) + 
      facet_wrap(vars(Scenarios),
                 nrow=12, ncol=3, scales="free_y", shrink=FALSE)

p10 <- ggplot() +
  geom_col(data = tePC, aes(x= TE, y=Count)) + 
  facet_wrap(vars(Scenarios),
             nrow=12, ncol=3, scales="free_y", shrink=FALSE)
p9

p10

2.2 Plots of finer population segments and model diagnostics

load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_longrun.rdata")
load("data/data.rdata")
library ('MCMCvis')
library ('coda')
library ('ggplot2')
library('reshape2')
library('bayestestR')
out <- lapply(post, as.mcmc)

# Identify chains with NAs that 
# failed to initialize
NAlist <- c()
for (i in 1:length(out)){
  NAlist[i] <- any (is.na(out[[i]][,1:286]) | out[[i]][,1:286]<0)
}
# Subset chains to those with good initial values
out <- out[!NAlist]
post2 <- post[!NAlist]
outp <- MCMCpstr(out, type="chains")

!NAlist
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
# default settings for plots 
plt  <- function(object, params,...) {
  MCMCplot(object=out, 
           params=params, 
           guide_axis=TRUE, 
           HPD=TRUE, ci=c(80, 95), horiz=FALSE, 
           #ylim=c(-10,10),
           ...)
  }

Plot model estimates of demographic rates. Life Stages are abbreviated as B = breeder, NB = nonbreeder, FY = first year. First-year abundance accounts for translocated birds.

# Abundance of females at Los Haitises
par(mfrow=c(5,2))
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NFY[",1:13, ", 1]"), 
    main="First-year (FY) minus hacked\n Los Haitises", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,1], 
     ylab="Counts", xlab="Year", type="b")

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NF[",1:13, ", 1]"), 
    main="Adult nonbreeder (NB)\n Los Haitises", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot.new()

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NB[",1:13, ", 1]"),
    main="Adult breeder (B)\n Los Haitises", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot.new()

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NAD[",1:13, ", 1]"), 
    main="Adult Breeders and Nonbreeders\n Los Haitises", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsAdults[,1], 
     ylab="Counts", xlab="Year",  type="b")

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("Ntot[",1:13, ", 1]"), 
    main="All stages\n Los Haitises", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,1]+datl$countsAdults[,1], 
     ylab="Counts", xlab="Year",  type="b")

# Abundance of females at Punta Cana
par(mfrow=c(5,2))
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NFY[",1:13, ", 2]"), 
    main="First-year (FY) plus hacked\n Punta Cana", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,2], 
     ylab="Counts", xlab="Year",  type="b")

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NF[",1:13, ", 2]"), 
    main="Adult nonbreeder (NB)\n Punta Cana", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot.new()

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NB[",1:13, ", 2]"),
    main="Adult breeder (B)\n Punta Cana", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot.new()

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NAD[",1:13, ", 2]"), 
    main="Adult Breeders and Nonbreeders\n Punta Cana", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsAdults[,2], 
     ylab="Counts", xlab="Year",  type="b")

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("Ntot[",1:13, ", 2]"), 
    main="All stages\n Punta Cana", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,2]+datl$countsAdults[,2], 
     ylab="Counts", xlab="Year",  type="b")

Population dynamics are determined by transitions, These transitions between stages are abbreviated as the starting life stage to the final life stage. For example a first-year recruiting to a breeder would be abbreviated as “FY to B”. I’ll list them here for convenience:

“FY to NB” is first-year to nonbreeder.

“NB to NB” is nonbreeder adult to nonbreeder adult.

“B to NB” is a breeding adult to a nonbreeder adult.

“FY to B” is first-year to breeder.

“NB to B” is nonbreeder adult to breeder adult.

“B to B” is breeder adult to breeder adult.

# Finer population segments
par(mfrow=c(4,2))
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 1, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nFirst-years fledged", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 2, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nFY to NB", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 3, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nNB to NB", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 4, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nB to NB", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 5, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nFY to B", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 6, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nNB to B",
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 7, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nB to B", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")


par(mfrow=c(4,2))

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 1, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nFY fledged",
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 2, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nFY to NB", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 3, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nNB to NB", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 4, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nB to NB", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 5, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nFY to B", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 6, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nNB to B",
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 7, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nB to B", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")

Other parameter estimates.

# I needed to abbreviate to save plot space
# FY=first-year, NB=Nonbreeder, B=Breeder
par(mfrow=c(1,2))
plt(object=out, 
    params=paste0("mus[",1:8, ", 1]"), 
    exact=TRUE, ISB=FALSE, 
    ylim=c(0,1),
    main="Overall means\n Los Haitises", 
    labels=c("FY survival", "NB survival", "B survival",
             "FY to B", "NB to B", "B to NB",
             "NB detection", "B detection")
    )

plt(object=out, 
    params=paste0("mus[",1:8, ", 2]"), 
    exact=TRUE, ISB=FALSE, 
    ylim=c(0,1),
    main="Overall means\n Punta Cana", 
    labels=c("FY survival", "NB survival", "B survival",
             "FY to B", "NB to B", "B to NB",
             "NB detection", "B detection"))

par(mfrow=c(1,1))
plt(object=out, 
    params="betas", 
    main= "Translocation effects",
    labels=c("FY survival", "NB survival", "B survival",
             "FY to B", "NB to B", "B to NB",
             "NB detection", "B detection"))

# Plot effort effects
plt(object=out, 
    params=paste0("deltas[",1:4,"]"),
    exact=TRUE, ISB=FALSE,
    main= "Effort effects", 
    labels=c("Number of Adults\nEffort", "Number of Adults\nEffort squared",
             "Number of FYs\nEffort", "Number of FYs\nEffort squared"),
    ylim=c(-0.5, 1))

plt(object=out, 
    params=paste0("deltas[",5:8,"]"),
    exact=TRUE, ISB=FALSE,
    main= "Effort effects continued", 
    labels=c("Nonbreeder detection\nEffort", "Nonbreeder detection\nEffort sq",
             "Breeder detection\nEffort", "Breeder detection\nEffort sq"),
    ylim=c(-5, 5))

# Fecundity
par(mfrow=c(1,1))
plt(object=out, 
    params=c("lmu.prod"), 
    labels= c("Productivity\n(log scale)\nLos Haitises",
              "Productivity\n(log scale)\nPunta Cana"))

f <- exp(outp$lmu.prod)

# Is fecundity at LHNP greater than PC
par(mfrow=c(1,1))
fdiff <- f[1,]-f[2,]
hist(fdiff, main="Productivity difference")
abline(v=0, lty=2)

# print probability of direction, similar to frequentist p-value
# so values <=0.025 and >=0.975
# Is the difference in fecundity >0 ?
mean(fdiff>0) 
## [1] 0.9988
# How many times greater is fecundity 
# at treated versus non-treated sites
# median(f.pred[1,]/f[1,])
# median(f.pred[2,]/f[2,])


# gamma = nest treatment effect on fecundity
par(mfrow=c(1,1))
plt(object=out, 
    params=c("gamma"), 
    main="Anti-Parasitic Fly\nTreatment Effects", ylim=c(0,3))

par(mfrow=c(1,1))
# sds <- paste0("sds[", 1:9, "]")
# plt(object=out, params=sds,
#     exact=TRUE, ISB=FALSE,
#     main="Temporal SDs (synchrony among sites)",
#     labels=c("FY survival", "NB survival", "B survival",
#              "FY to B", "NB to B", "B to NB",
#              "NB detection", "B detection",
#              "Productivity"))
sds2 <- paste0("sds2[", 1:9, "]")
plt(object=out, params=sds2,
    exact=TRUE, ISB=FALSE,
    main="Site-temporal SDs", 
    labels=c("FY survival", "NB survival", "B survival",
             "FY to B", "NB to B", "B to NB",
             "NB detection", "B detection",
             "Productivity"))

# Correlations among vital rates
# Plot is messy with only a few strong correlations
ind <- 1
Rs <- R2s <- c()
for (i in 1:(nrow(outp$R)-1)){
  for (j in (i+1):nrow(outp$R)){
  Rs[ind] <- paste0("R[",i,", ", j, "]")
  R2s[ind] <- paste0("R2[",i,", ", j, "]")
  ind <- ind+1
  }}
# par(mfrow=c(2,1))
# plt(object=out, params=Rs[1:18], exact=TRUE, ISB=FALSE,
#     main="Correlations btw demographic rates\n over time (synchrony)",
#     xlab = "Rhos", guide_lines=TRUE)
# plt(object=out, params=Rs[19:36], exact=TRUE, ISB=FALSE,
#     main="Correlations btw demographic rates\n over time (synchrony), continued...",
#     xlab = "Rhos", guide_lines=TRUE)
par(mfrow=c(2,1))
plt(object=out, params=R2s[1:18], exact=TRUE, ISB=FALSE, 
    main="Correlations btw demographic rates\n over time and sites",
    xlab = "Rhos", guide_lines=TRUE)
plt(object=out, params=R2s[19:36], exact=TRUE, ISB=FALSE, 
    main="Correlations btw demographic rates\n over time and sites, continued ...",
    xlab = "Rhos", guide_lines=TRUE)

# Annual averages for integration into the population model
par(mfrow=c(1,1))
labs <- c(paste0("LH ",2011:2023), paste0("PC ",2011:2023))
plt(object=out, params="mn.phiFY", ylim=c(0,1),
    main="First-year survival", labels = labs,
    xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.phiA", ylim=c(0,1),
    main="Adult nonbreeder", labels = labs,
    xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.phiB", ylim=c(0,1),
    main="Breeder", labels = labs,
    xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.psiFYB", ylim=c(0,1),
    main="First-year to breeder", labels = labs,
    xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.psiAB", ylim=c(0,1),
    main="Adult nonbreeder to breeder", labels = labs,
    xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.psiBA", ylim=c(0,1),
    main="Adult breeder to nonbreeder", labels = labs,
    xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.pA", ylim=c(0,1),
    main="Nonbreeder", labels = labs,
    xlab = "Year", ylab= "Detection")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.pB", ylim=c(0,1),
    main="Breeder", labels = labs,
    xlab = "Year", ylab= "Detection")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.prod",
    main="", labels=labs,
    xlab = "Year", ylab= "Productivity")
abline(v=13.5, lwd=2)

3. Model diagnostics

3.1 Check Goodness-of-fit

Goodness-of-fit tests provide evidence that statistical distributions adequately describe the data. Here we test fit for brood size and counts. A Bayesian p-value nearest to 0.5 suggests good fitting statistical distributions, while values near 1 or 0 suggest poor fit.

# Goodness of fit check
fit.check <- function(out, ratio=FALSE, 
                      name.rep="f.dmape.rep", 
                      name.obs="f.dmape.obs",
                      jit=100,
                      ind=1,
                      lab=""){
  par(mfrow=c(1,1))
  # plot mean absolute percentage error
  samps <- MCMCpstr(out, "all", type="chains")
  rep <- samps[name.rep][[1]][ind,]
  obs <- samps[name.obs][[1]][ind,]
  mx <- max(c(rep, obs))
  mn <- min(c(rep, obs))
  plot(jitter(obs, amount=jit), 
       jitter(rep, amount=jit),
       main=paste0("Mean absolute percentage error\n",lab),
       ylab="Discrepancy replicate values",
       xlab="Discrepancy observed values", 
       xlim=c(mn,mx), ylim=c(mn,mx), 
       pch=16, cex=0.5, col="gray10")
  curve(1*x, from=mn, to=mx, add=T, lty=2, lwd=2, col="blue")
  bp1 <- round(mean(rep > obs),2)
  loc <- ifelse(bp1 < 0.5, "topleft", "bottomright")
  legend(loc, legend=bquote(p[B]==.(bp1)), bty="n", cex=2)
  
  if (ratio==TRUE){
    t.rep <- samps["tvm.rep"][[1]][ind,]
    t.obs <- samps["tvm.obs"][[1]][ind,]
    # plot variance/mean ratio
    hist(t.rep, nclass=50,
         xlab="variance/mean ", main=NA, axes=FALSE)
    abline(v=t.obs, col="red")
    axis(1); axis(2)
  }
  return(list('Bayesian p-value'=bp1))
}

# Counts of adults (nonbreeders+breeders)
fit.check(out, ratio=F,
          name.rep="dmape.rep", 
          name.obs="dmape.obs",
          ind=1,
          lab="Adults(Breeder+Nonbreeder)- Poisson", jit=300)

## $`Bayesian p-value`
## [1] 0.66
# Counts of first-year fledglings, ind=2
# Poisson failed fit test bp=0
# So we assigned negative binomial only for Punta Cana.
# Los Haitises excluded from neg bin because no zeroes in counts. 
fit.check(out, ratio=F,
          name.rep="dmape.rep", 
          name.obs="dmape.obs",
          ind=2,
          lab="First-year counts\nNeg binomial-Poisson", jit=300)

## $`Bayesian p-value`
## [1] 0.57
# Productivity
fit.check(out, ratio=F,
          name.rep="f.dmape.rep", 
          name.obs="f.dmape.obs",
          ind=1,
          lab="Productivity-Neg binomial", jit=300)

## $`Bayesian p-value`
## [1] 0.83

3.2 Examine Traceplots, Compare Posteriors with Priors

Traceplots provide evidence that models adequately converged.

MCMCtrace(post2, pdf=FALSE, params= "mus", Rhat=TRUE, priors=runif(20000, 0, 1), post_zm=FALSE)

MCMCtrace(post2, pdf=FALSE, params= "betas", Rhat=TRUE, priors=runif(20000, -20, 20), post_zm=FALSE)

MCMCtrace(post2, pdf=FALSE, params= "deltas", Rhat=TRUE, priors=runif(20000, -20, 20), post_zm=FALSE)

MCMCtrace(post2, pdf=FALSE, params= "gamma", Rhat=TRUE, priors=runif(20000, -20, 20), post_zm=FALSE)

MCMCtrace(post2, pdf=FALSE, params= "sds2", Rhat=TRUE, priors=rexp(20000, 1), post_zm=FALSE)

MCMCtrace(post2, pdf=FALSE, params= "R2")